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O ABSTRACT 

(N 

Mh We investigate convergence of Lagrangian Perturbation Theory (LPT) by analysing 

r~^ the model problem of a spherical homogeneous top-hat in an Einstein-deSitter back- 

ground cosmology. We derive the formal structure of the LPT series expansion, working 
\l to arbitrary order in the initial perturbation amplitude. The factors that regulate LPT 

convergence are identified by studying the exact, analytic solution expanded accord- 
I— I ing to this formal structure. The key methodology is to complexify the exact solution, 

^^ demonstrate that it is analytic and apply well-known convergence criteria for power 

r\) series expansions of analytic functions. The "radius of convergence" and the "time of 

validity" for the LPT expansion are of great practical interest. The former describes 
the range of initial perturbation amplitudes which converge over some fixed, future 
time interval. The latter describes the extent in time for convergence of a given ini- 
O tial amplitude. We determine the radius of convergence and time of validity for a full 

-J-J sampling of initial density and velocity perturbations. 

j^ This analysis fully explains the previously reported observation that LPT fails to 

'"^ predict the evolution of an underdense, open region beyond a certain time. It also 

^vj implies the existence of other examples, including overdense, closed regions, for which 

^ LPT predictions should also fail. We show that this is indeed the case by numerically 

t^^ computing the LPT expansion in these problematic cases. 

The formal limitations to the validity of LPT expansion are considerably more 

complicated than simply the first occurrence of orbit crossings as is often assumed. 

Evolution to a future time generically requires re-expanding the solution in overlapping 

JLj domains that ultimately link the initial and final times, each domain subject to its 

^— s own time of validity criterion. We demonstrate that it is possible to handle all the 

,—1 problematic cases by taking multiple steps (LPT re-expansion). 

L| A relatively small number (~ 10) of re-expansion steps suffices to satisfy the time 

of validity constraints for calculating the evolution of a non-collapsed, recombination- 
era perturbation up to the current epoch. If it were possible to work to infinite La- 
grangian order then the result would be exact. Instead, a finite expansion has finite 
errors. We characterise how the leading order numerical error for a solution generated 
by LPT re-expansion varies with the choice of Lagrangian order and of time step size. 
Convergence occurs when the Lagrangian order increases and/or the time step size 
decreases in a simple, well-defined manner. We develop a recipe for time step control 
for LPT re-expansion based on these results. 

Key words: cosmology: theory ~ large-scale structure of Universe. 



* E-mail: sirLn27@cornell.edu 

t E-mail: cliornoff@astro.cornoll.edu 

© RAS 



2 Sharvari Nadkarni- Ghosh and David F. Chernoff 

1 INTRODUCTION 

Understanding the non-linear growth of structure in an expanding universe has been an active area of research for nearly 
four decades. Simulations have been instrumental in illustrating exactly what happens to an initial power spectrum of small 
fluctuations but analytic methods remain essential for elucidating the physical basis of the numerical results. Perturbation 
theory, in particular, is an invaluable tool for achieving a sophisticated understanding. 

The Eulerian and Lagrangian frameworks are the two principal modes of description of a fluid. The fundamental dependent 
variables in the Eulerian treatment are the density p(x, t) and velocity v(x, t) expressed as functions of the grid coordinates 
X and time t, the independent variables. In perturbation theory the dependent functions are expanded in powers of a small 
parameter. For cosmology that parameter typically encodes a characteristic small spatial variation of density and/or velocity 
with respect to a homogeneous cosmology at the initial time. As a practical matter, the first-order perturbation theory 
becomes inaccurate when the perturbation grows to order unity. Subsequently one must work to higher order to handle the 



development of non-linearity (see Bernardeau et al. 112002 for a review) or adopt an alternative method of expansion. 



In the Lagrangian framework, the fundamental dependent variable is the physical position of a fluid element or particle 
(terms used interchangeably here). The independent variables are a set of labels X, each of which follows a fluid element, and 
the time. Usually X is taken as the position of the element at some initial time but other choices are possible. In any case, 
the physical position and velocity of a fluid element are r = r(X,t) and r(X,i), respectively. Knowledge of the motion of 
each fluid element permits the full reconstruction of the Eulerian density and velocity fields. In cosmological applications of 
Lagrangian perturbation theory (LPT), just like Eulerian perturbation theory, the dependent variables are expanded in terms 
of initial deviations with respect to a homogeneous background. The crucial difference is that the basis for the expansion 
is the variation in the initial position and position-derivative not the variation in the initial fiuid density and velocity. The 
Eulerian density and velocity may be reconstructed from knowledge of the Lagrangian position using exact non-perturbative 
definitions. A linear approximation to the displacement field results in a non-linear expression for the density contrast. The 
Lagrangian description is well-suited to smooth, well-ordered initial conditions; a single fiuid treatment breaks down once 
particle crossings begin, caustics form and the density formally diverges. 



First-order LPT was originally introduced by Zel'Dovich ( 1970 1 to study the formation of non-linear structure in cos 



mology. In his treatment the initial density field was taken to be linearly proportional to the initial displacement field (the 
"Zeldovich approximation"). These results were extended by many authors ([Moutarde et al.|1991 Buchert|1992[ Bouchet et al 



Buchcrt & Ehlcrs 1993; Buchert][T994] [Munshi et aL][T994l |Catelan|[T995l |Buchert|[T995l [Bouchet et"aI][T995l [Bouchet 
.Ehlers fc Buchcrt,|1997,). The work pioneered by Bouchet focused on Zeldovich initial conditions and established the 



1992 
T996, 



link between LPT variables and statistical observables. The work by Buchert as well as the paper by Ehlers & Buchert (19971 



formalised the structure of the Newtonian perturbative series for arbitrary initial conditions. A general relativistic version of 



the Zeldovich approximation was developed by Kasai ( 1995 1 and other relativistic descriptions of the fiuid in its rest frame 
were investigated by Matarrese & Terranova (19961 and Matarrese, Pantano & Saez (1993 [1994[). LPT has been used for 



many applications including, recently, the construction of non-linear halo mass functions by Monaco| ( |1997| ) and [Scoccimarro] 
& Sheth[p002l). 



Not much has been written about the convergence of LPT although LPT expansions are routinely employed. [Sahni fc] 
[Shandann] ( |1996| pointed out that the formal series solution for the simplest problem, the spherical top-hat, did not converge 
for the evolution of homogeneous voids. Figure [l] illustrates the conundrum that the LPT approximations diverge from the 
exact solution in a manner that worsens as the order of the approximation increases. The details will be described in the next 
section. 

This paper explores LPT convergence for the spherical top-hat and identifies the root cause for the lack of convergence. 
The analysis naturally suggests a means of extending the range of validity of LPT. This generalisation of LPT guarantees 
convergence to the exact solution of the model problem at all times prior to the occurrence of the first caustic. 

[Tatekawa] ( |2007| ) attempted to treat the divergence by applying the Shanks transformation to the LPT series. Although 
non-linear transformations can sum a divergent series, the correct answer is not guaranteed; comparison of several different 
methods is usually necessary to yield trustworthy results. Other approaches include the Shifted-Time-Approximation (STA) 



and Frozen-Time- Approximation (FTA) which have been investigated by Karakatsanis, Buchert & Melott (19971. These 



schemes modify lower order terms to mimic the behavior of higher order terms and/or extend the range of applicability in 
time. None of these techniques are considered here. 

The organisation follows: ^sketches the model problem, the evolution of a uniform sphere in a background homogeneous 
Einstein-deSitter cosmology. The LPT equations, the structure of the formal series and the term-by-term solution are outlined. 
33] discusses the complexification of the LPT solution and convergence of the series. This section introduces the "radius of 
convergence" and the "time of validity" for LPT. 34] outlines the real and complex forms of the parametric solution and sets 
forth the equations that must be solved to locate the poles which govern the convergence. i|5]presents numerical results for the 
time of validity and radius of convergence for a full range of possible initial conditions for the top-hat. The notion of mirror 
model symmetry is introduced and used to explain a connection in the convergence for open and closed models. !|6] shows that 
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Figure 1. The time-dependent scale factor b of an initial spherical top-hat perturbation is plotted as a function of the background scale 
factor a. The perturbation is a pure growing mode, i.e. the density and velocity perturbations vanish at i = 0. The black dotted line is 
the exact solution. The smooth blue lines are the LPT results obtained by working successively to higher and higher order. Series with 
even (odd) final order lie below (above) the exact solution. Roughly speaking, LPT converges only for a < 0.2. Beyond that point the 
higher order approximations deviate from the exact solution more than lower order ones. 

the time of validity may be extended by re-expanding the solution in overlapping domains that ultimately link the initial and 
final times, each domain subject to an individual time of validity criterion. The feasibility of this method is demonstrated in 
some examples, ^summarises the work. 



2 THE MODEL PROBLEM AND FORMAL SERIES SOLUTION 

This section describes the governing equations, the initial physical conditions, the formal structure of the LPT series solution 
and the order-by-order solution. 

2.1 Newtonian treatment 

Consider evolution on sub-horizon scales after recombination in a matter-dominated universe. A Newtonian treatment of 
gravity based on solving Poisson's equation for the scalar potential and on evaluating the force in terms of the gradient of 
the potential gives an excellent approximation for non-relativistic dynamics. When there are no significant additional forces 
on the fluid element (e.g. pressure forces) then it is straightforward to eliminate the gradient of the potential in favour of i", 
the acceleration. The governing equations are 






-AnGp(TC,t) 







(1) 
(2) 



where p(x, t) is the background plus perturbation density, G is Newton's gravitational constant and Va; is the Eulerian gradient 
operator. In the Lagrangian treatment, the independent variables are transformed (x, t) — >■ (X, t) and the particle position 
r = r(X, t) adopted as the fundamental dependent quantity. For clarity note that x refers to a fbced Eulerian grid not a 
comoving coordinate. 

2.2 Spherical top- hat 

The starting physical configuration is a compensated spherical perturbation in a homogeneous background cosmology. The 
perturbation encompasses a constant density sphere about the centre of symmetry and a compensating spherical shell. The 
shell that surrounds the sphere may include vacuum regions plus regions of varying density. Unperturbed background extends 
beyond the outer edge of the shell. Physical distances are measured with respect to the centre of symmetry. At initial time to 
the background and the innermost perturbed spherical region (hereafter, "the sphere") have Hubble constants Ho and Hpo, 
and densities po and ppo, respectively. Let rbfi (rpfi) be the physical distance from the centre of symmetry to the inner edge of 
the background (to the outer edge of the sphere) at the initial time. Let ao, bo be the initial scale factors for the background 
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and the sphere respectively. Two sets of Lagrangian coordinates Y = rbfi/ao and X = rp^/bo are defined. A gauge choice sets 
ao = bo- Appendix IA] provides a figure and gives a somewhat more detailed chain of reasoning that clarifies the construction 
of the physical and Lagrangian coordinate systems. The initial perturbation is characterised by the independent parameters 

S = ^-1 
Po 

5. = t^-1- (3) 

Finally, assume that the background cosmology is critical Slo = 1. The perturbed sphere has 

"^o=(TTMl- (4) 

The physical problem of interest here is the future evolution of an arbitrary initial state unconstrained by the past history. In 
general, the background and the perturbation can have different big bang times. Initial conditions with equal big bang times 
will be analysed as a special case of interest and imply an additional relationship between 5 and 5^ . 

While the previous paragraphs summarise the set up, they eschew the complications in modelling an inhomogeneous 
system in terms of separate inner and outer homogeneous universes. For example, matter motions within the perturbed inner 
region may overtake the outer homogeneous region so that there are problem-specific limits on how long solutions for the 
scale factors a{t) and b{t) remain valid. The appendix shows that there exist inhomogeneous initial configurations for which 
the limitations arising from the convergence of the LPT series are completely independent of the limitations associated with 
collisions or crossings of inner and outer matter-filled regions. A basic premise of this paper is that it is useful to explore the 
limitations of the LPT series independent of the additional complications that inhomogeneity entails. 

2.3 Equation governing scale factors 

During the time that the spherical perturbation evolves as an independent homogeneous universe it may be fully described 
in terms of the motion of its outer edge Vp . Write 

rp{t) = b{t)X (5) 

where b{t) is the scale factor and X is the Lagrangian coordinate of the edge. The initial matter density of the homogeneous 
sphere p{X, to) = ppo = po{l + S). The physical density of the perturbation at time t is 

^rv -f\ piX,to)J(X,to) 

P^^''^ = — JixJ) — (^) 

where the Jacobian of the transformation relating the Lagrangian and physical spaces is 

J{X,t)=det(^^y (7) 

Since eq. ([5| implies J{X,t) — b{t)^ and the choice ao = &o implies J{X,U)) — qq the perturbation matter density at later 
times is 

_ po{l + 5)al 

p^^^> - b{tr ■ ^' 

Substituting for pp and rp in eq. ([l| gives 

b l Hial{l + 6) 

6^2 &3 



(9) 



with initial conditions fe(to) = lo and 6(fo) ~ ao(l + 5„). The curl of the acceleration (i.e. eq. ([2|) vanishes by spherical 
symmetry. The corresponding equation for the background scale factor is 

a ^ 1 HqUq 

a 2 a3 ^ ' 

with initial conditions a{to) — oq and a(to) = ao = aoHo. The solution for b{t) will be expressed in terms of its deviations 
from a{t). 

In summary, the physical setup is an fio = 1 background model and a compensated spherical top-hat (over- or underdense). 
The properties of interest are the relative scale factors a(t)/ao and b{t)/ao (the choice of ao is arbitrary and bo = ao). The 
evolution of the relative scale factors is fully specified by _ffo, ^po and VLpo at time to- The perturbed physical quantities, Hpo 
and Slpo, may be equivalently specified by a choice of 5 and 5„. Appendix lA] contains a systematic description and enumerates 
degrees of freedom, parameters, constraints, etc. 
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Figure 2. Phase diagram of density and velocity perturbations (<5, <5„). Physical initial conditions require —1 < 5 < oo and — oo < 5„ < oo. 
The left panel highlights the qualitatively different initial conditions. The shaded (unshaded) region corresponds to closed (open) model 
with negative (positive) total energy. For small A, models with 9^ < 8 < 9^ are closed. Initially expanding and contracting models are 
separated by the dashed horizontal line {5y = —1). The right panel shows the evolution of S and 5y. The solid blue line corresponds to 
the "Zeldovich" condition i.e no perturbation at t = 0. The points (— 1, — 1), (— 1, 0.5) and (0,0) are unstable, stable and saddle fixed 
points of the phase space flow. The flow lines (indicated by the blue vectors) converge along the Zeldovich curve either to the stable fixed 
point at (—1,0.5) or move parallel to the Zeldovich curve to a future density singularity. Further discussion follows in i]5.4|and j]6.2[ 



2.4 Perturbations in phase space 



The initial density and velocity perturbations are taken to be of the same order in the formalism developed by [Buchert] 
(1992) |1994[), [Buchert fc Ehlers| (|1993|) and|Ehlers fc Buchert] (|1997[). We assume the same ordering here. Write the initial 



perturbation {S, 5„) in terms of magnitude A and angle 9 



A = v/52 + 52 

so that 

5 = A cos 61 
5v = Asin6'. 



(11) 



(12) 
(13) 



To map physical perturbations {5,Sv) in a unique manner to (A, 6*) adopt the ranges A Js and — tt < 6 ^ it. Figure [2] (left 
panel) shows the phase space of initial perturbations. Since density is non-negative the regime of physical interest is S ^ —1. 
Open (closed) models with positive (negative) total energy are the regions that are unshaded (shaded). Initially expanding 
models, 1 + 5y > 0, lie above the horizontal dashed line. The right panel of figure [2] summarises the overall evolution of the 
system. The initial choice of S and 5„ dictates the trajectory in the plane. Cosmologically relevant initial conditions generally 
assume there to be no perturbation at f = 0. We adopt the name "Zeldovich" initial conditions for models that satisfy this 
condition. This establishes a specific relation between 5 and Sv which is indicated by the sold blue line. The exact mathematical 



relationship is given in f5.4 Starting from a general initial point (5, 5„), the system as it evolves traces out a curve in phase 
space indicated by the blue arrows. There are three fixed points visible. The origin {S,5v) = (0,0), which corresponds to a 
unperturbed background model, is a saddle point. The vacuum static model at point (—1,-1) is a unstable node and the 
vacuum, expanding model at (—1,0.5) is a degenerate attracting node. Far to the right and below the dashed line the models 
collapse to a future singularity. The phase portrait illustrates that the trajectories either converge to the vacuum, expanding 
model or to the singular, collapsing model. The equations that govern the flow and further relevance of the Zeldovich solution 
is discussed in 96.21 and 35.41 
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2.5 Generating the Lagrangian series solution 

The scale factor is formally expanded 

where fe'"' denotes an n-th order term. The initial conditions are 

fo(io) = a{to) 

b{io) = ao{l + Sy)^ho{l + Asme). 

Substitute the expansion for b{t) into eq. (IM, equate orders of A to give at zeroth order 



h(0) 



1 Hpao 

2 6(0)2 



= 



(14) 



(15) 
(16) 



(17) 



which is identical in form to eq. 1 10 1 for the unperturbed background scale factor. The initial conditions at zeroth order: 

fo'°'(io) = ao (18) 

fo™(io) = ao. (19) 

The equation and initial conditions for 6''^'(f) simply reproduce the background scale factor evolution b^'^'{t) = a{t). Without 
loss of generality assume that the background model has big bang time i = so that 

At first order 



hW 



H^alb'-^'' 1 Rial cose 



and, in general, 
Hialb^'^^ 



l(") 



-(") 



(21) 



(22) 



where S^"' depends upon lower order approximations (6'° , b^^' . . . ft'" ^') as well as 9. The first few are: 



^(2) ^ _1:(M^ [;,(!) {35(1) „ 2a cos e}] 



2 a* 

(3) _ 1 HqOq •" 

2 flS 

(4) _ 1 ^0 fflp ■" 



20^6'^^ cos ( 



2 aS 



,W|_4(6(^')' + 6afe(2)+3a6«cos0} 
(6(^))' {5 (b(^')' - 12a6(2) - 4a6(^' cos^} + 
Ba^b'^' {&(^) + 6<2' cose} + Sa^ (b'^')' - 2a^b(^) cos( 



These terms can be easily generated by symbolic manipulation software. The initial conditions are 

b^'\to) = 

b (to) ~ fflo sin0 

and for n > 1 

b^"Hto) = 

b^"\to) = 0. 



(23) 
(24) 

(25) 



(26) 
(27) 

(28) 
(29) 



The ordinary differential equations for b^"' may be solved order- by-order. 

To summarise, the structure of the hierarchy and the simplicity of the initial conditions allows the evaluation of the 
solution at any given order in terms of the solutions with lower order. This yields a formal expansion for the scale factor of 
the sphere 



b = ^b'"'(i)A" 

n=0 



(30) 
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Extending the domain of validity of the Lagrangian approximation 7 

which encapsulates the Lagrangian perturbation treatment. The right hand size explicitly depends upon the size of the 
perturbation and time and implicitly upon ao, Hq, and 8. This hierarchy of equations is identical to that generated by the 
full formalism developed by Buchert and collaborators when it is applied to the top-hat problem. The convergence properties 
in time and in A are distinct; a simple illustrative example of this phenomenon is presented in Appendix [B] 



3 CONVERGENCE PROPERTIES OF THE LPT SERIES SOLUTION 

The series solution outlined in the previous section does not converge at all times. Figure IT] is a practical demonstration of 
this non-convergence for the case of an expanding void. An understanding of the convergence of the LPT series is achieved 
by extending the domain of the expansion variable A from the real positive axis to the complex plane. 

3.1 Complexification 

The differential eq. ([9l and initial conditions for the physical system are 

■■ _ l iJgag(l-HAcose) 

^' ' 2 b{ty 

b{to) = ao 

fe(io) = do(l + Asin6») (31) 

where t, b{t), A and all zero-subscripted quantities are real. This set may be extended by allowing A and b to become complex 
quantities, denoted hereafter, A and b, while the rest of the variables remain real. The complex set is 

^^) _ IHialil + ^cose) 

b(io) 

b(io) = ao(l-h Asine*). (32) 



2 


b(t)2 


ao 




ao(l + A 


sine*). 



The theory of differential equations (for example, Chicone 2006 1 guarantees that the solution to a real initial value problem 



is unique and smooth in the initial conditions and parameters of the equation and can be extended in time as long as there are 
no singularities in the differential equation (hereafter, the maximum extension of the solution). First, note that each complex 



quantity in eq. (32 1 may be represented by a real pair, i.e. h = u + iv hy pair {u,v} = {5Rb, SJb} and A = a; -I- ij/ by pair 
{x, y} — {5RA, SJA}. The basic theory implies continuity and smoothness of solution u and v with respect to initial conditions 
and parameters x and y. Second, observe that the Cauchy-Riemann conditions Ux = Vy and Uy = —Vx are preserved by the 
form of the ordinary differential equation. Since the initial conditions and parameter dependence are holomorphic functions 
of A it follows that b(i. A) is a holomorphic function of A at times t within the maximum extension of the solution. 

Inspection shows that the differential equation is singular only at b = 0. For a particular value of A = A', the solution 
to the initial value problem can be extended to a maximum time t^x such that b(A', tmx) = or to infinity. The existence of 
a finite tmx signals that a pole in the complex analytic function b(A,i) forms at A = A' and t = tmx- For times t such that 
to ^ t < tmx, the solution 6(A,i) is analytic in a small neighbourhood around the point A'. Of course, there may be poles 
elsewhere in the complex A plane. 

The relationship between the original, real-valued physical problem and the complexified system is the following. In the 
original problem A is a real, positive quantity at to- LPT is a power series expansion in A about the origin (the point A = 0). 
LPT's convergence at any time i can be understood by study of the complexified system. Consider the complex disk D centred 
on the origin and defined by |A| < A. At io each point in D determines a trajectory b(A,i) for the complexified system 
extending to infinity or limited to finite time t — tmi(A) because of the occurrence of a pole. The time of validity is defined 
as r(A) — mini7 tmx, i.e. the minimum tmx over the disk. Since there are no poles in D at to the time of validity is the span 
of time when D remains clear of any singularities. If a function of a complex variable is analytic throughout an open disk 
centred around a given point in the complex plane then the series expansion of the function around that point is convergent 



(Brown & Churchill 19961. The LPT expansion for the original problem converges for times less than the time of validity 
because the complex extension b(A,f) is analytic throughout V for t < T(A). If Ai < A2 then, in an obvious notation, the 
disks are nested I?(Ai) C 'D{A2) and the times of validity are ordered r(Ai) ^ T(A2). 

This idea is shown in figure IS] No singularities are present for the initial conditions at to; at ti a singularity is present 
outside the disk but it does not prevent the convergence of the LPT expansion with A equal to the disk radius shown; at t2 
a singularity is present in the disk or on its boundary and it may interfere with convergence. 

A distinct but related concept is the maximum amplitude perturbation for which the LPT expansion converges at the 
initial time and at all intermediate times up to a given time. The radius of convergence i?A(t) is the maximum disk radius A 
for which t > T{A). Because the disks are nested if ti < t2 then i?A(ti) ^ RA{t2)- 
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t = t„ 



t = t, 



t = t. 



Figure 3. This figure is a schematic illustration of how the time of validity is determined. The initial conditions imply a specific, real 
A at time to- The LPT series is an expansion about A = 0, convergent until a pole appears at some later time within the disk of radius 
A (shown in cyan) in the complex A plane. Typically, the pole's position forms a curve (blue dashed) in the three dimensional space 
(5R[A], S[A], t). The black dots mark the pole at times t\ and t2. At t\ the pole does not interfere with the convergence of the LPT series; 
at t2 it does. The time of validity may be determined by a pole that appears within the disk without moving through the boundary (not 
illustrated). 

The time of validity and the radius of convergence are inverse functions of each other. If the initial perturbation is 
specified, i.e. A is fixed, and the question to be answered is "how far into the future does LPT work?" then the time of 
validity gives the answer. However, if the question is "how big an initial perturbation will be properly approximated by LPT 
over a given time interval?" then the radius of convergence provides the answer. 

Finally, note that one can trivially extend this formalism to deal with time intervals in the past. 

3.2 Calculating radius of convergence and time of validity 

The following recipe shows how to calculate the radius of convergence R^{t) and the time of validity T(A) efficiently. Fix ao. 
Ho, to and 9; these are all real constants set by the initial conditions. Assume that it is possible to find b(A,f) for complex 



A and real t by solving eq. (321. There exist explicit expressions for b as will be shown later. 



Start with t = to and R^{t) — oo. The iteration below maps out RA{t) by making small increments in time St. 



Store old time ip 



= i, choose increment St and form new time of interest t = tprevious + St. 



• Locate all the A which solve b(A,f) = 0. The roots correspond to poles in the complex function. Find the root closest 
to the origin and denote its distance as |Anear|. 

• The radius of convergence is -RA(t) ~ min(|A„ear|, flA(tpreDioiis)). 

• Continue. 

Since Ra is decreasing, the inversion to form T{A) is straightforward. Figure W] shows a schematic cartoon of the con- 
struction process. 
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Figure 4. A schematic illustration of the radius of convergence and the time of validity. The left panel shows the location of poles in 
the complex A plane at times ti and t2, denoted by orange squares and green dots, respectively. At a fixed time, the pole nearest the 
origin determines the disk (black circle) within which a series expansion about the origin converges. The right panel shows |A| for ti 
and i2- The black line is R^{t), the minimum |A| calculated for a continuous range of times (where to, the initial time, lies far to the 
left). The arrows show how the time of validity is inferred for a given A. 

4 EXPLICIT SOLUTIONS 

The usual parametric representation provides an efficient method to construct an explicit complex representation for b(A,t). 



4.1 Real (physical) solutions 



The original system eq. (31 1 depends upon ao, Ho, 9 and A. The assumed Einstein-deSitter background has ao > and 



do > 0; as defined, the perturbation amplitude A Js and the relative density and velocity components are determined by 
phase angle 6 with — tt < 9 ^ n. The quantity (1 + Acos6') is proportional to total density and must be non-negative. The 
sign of bo is the sign of 1 + Asin^ and encodes expanding and contracting initial conditions. 
Briefiy reviewing the usual physical solution, the integrated form is 

(1 + Asin^ 



(1 + Acos( 
6 



+ 



^-(1 + Acos^ 



The combination 



E{A,9) = (1 + Asin^)^- (1 + Acos^ 



(33) 



(34) 



is proportional to the total energy of the system. If £ > the model is open and if _B < it is closed and will re-collapse 
eventually. Figure [2] shows the parabola E = which separates open and closed regions. For infinitesimal A the line of division 
has slope tan 6* = 1/2. Models with 9 G [9^,9+] = [-7r-|-tan^^(l/2), tan~^(l/2)] = [-2.68,0.46] are closed while those outside 
this range are open. 

There are four types of initial conditions (positive and negative E, positive and negative bo) and four types of solutions, 
shown schematically in figure [5] The solutions have well-known parametric forms involving trigonometric functions of angle rj 
or iri (see Appendix[c|. The convention adopted here is that the singularity nearest the initial time io coincides with r; — and 



is denoted t 
and to is ta 



itl 



hang \ bang 



for initially expanding (contracting) solutions (see figure 5|. The time interval between the singularity 



ge 



|io-tr.„„1^0 
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Figure 5. Scale factor as a function of time. The initial conditions (60 = no = 1 and varying 60) are given at time to (dashed blue line). 
The left (right) panel illustrates initially expanding (contracting) models, t, corresponds to »; = 0; ti-di to rj = 1-k. For expanding 

solutions tage = ^0 ~ tjT is the time interval since the initial singularity and t^all is the future singularity for closed models. For 



contracting solutions tags 



bang 



to is the time until the final singularity and tcoU is the past singularity for closed models. 



The parametric solution for the models can be written as 

./ A a\ "-o (1 + Acos6l) 
b(r},A,0) — — ^7 — -— ^(l — cosri) 



t{7j,A,e) 



, , 1 (1 + Acose) , . , 

^0± 1 TTTT-. ^-. ZZTT^iV - Smri) - tage[A, 



2Ho [-E{A,e)] 



3/2 ' 



(35) 



The plus and minus signs give the solution for initially expanding and initially contracting models respectively. Parameter r] 
is purely real for closed solutions and purely imaginary for open solutions. The distance to the nearest singularity is 



db 



■y=i 



dy 



The second equality uses eq. (|33[) and the substitution y = b/ao. 



(36) 



4.2 Complex extension 

To extend the above parametric solution to the complex plane, one might guess the substitution A — >■ Ae"^ where —tv < 4> ^ tt 



in eq. (351 and eq. (36 1. The physical limit is <j!> = 0. However, this leads to two problems. First, the integral for tage can have 



multiple extensions that agree for physical (f> = but differ elsewhere including the negative real axis. This is tied to the fact 
that the operations of integration and substitution A — > Ae"^ do not commute because of the presence of the square root in 
the expression for tage- A second related problem is the presence of multiple square roots in the parametric form for t. These 
give rise to discontinuities along branch cuts such that one parametric form need not be valid for the entire range of (j), but 
instead the solution may switch between different forms. Directly extending the parametric solution is cumbersome. 



However, the original differential eq. (321 is manifestly single- valued. The equation can be integrated forward or backward 



numerically to obtain the correct solution for complex A. One can then match the numerical solution to the above parametric 
forms to select the correct branch cuts. This procedure was implemented to obtain the form for all A and 9. The main result 
is that the solution space for all 6 and A is completely spanned by complex extensions of the two real parametric forms which 
describe initially expanding and contracting solutions. The expressions for tage and details are given in Appendix |C3[ 

The traditional textbook treatment relating physical cosmological models with real fi > 1 and SI < 1 typically invokes 
a discrete transformation 77 — )■ ir; in the parametric forms and one verifies that this exchanges closed and open solutions. 
However, starting from the second order differential equation it is straightforward to use the same type of reasoning as above 
to construct an explicit analytic continuation from one physical regime to the other. 
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In addition, note that the differential equation and its solution remain unchanged under the simultaneous transformations 
A — > —A and 9 ^^ 6 + n. Every complex solution with — vr < ^ ^ can be mapped to a complex solution with < 6 ^ tt 
and vice- versa. For determining the radius of convergence and the time of validity the whole disk of radius | A| is searched for 
poles so it suffices to consider a restricted range of 6 to handle all physical initial conditions. 

4.3 Poles 

The condition b = signals the presence of a pole. Inspection of the parametric form shows that this condition can occur 
only when 77 = or r; = 27r. The corresponding time 

t(A,.) = /^o±(^^^±^^-W(A)) (, = 2.) ^3^^ 

t i0=Ftage(A) (»7 = 0) 

is immediately inferred. Since the independent variable t is real the transcendental equation 

9t(A,6l) = (38) 

must be solved. It is straightforward to scan the complex A plane and calculate t to locate solutions. Each solution gives a 
root of b = and also implies the existence of a pole at the corresponding A. Note that relying upon the parametric solutions 
is a far more efficient method for finding the poles than integrating the complex differential equations numerically. We have 
verified that both methods produce the same results. 

In practice, we fix 9, scan a large area of the complex A plane, locate all purely real t and save the {A, t} pairs. These 
are used to create a scatter plot of | A| as a function of time (hereafter the "root plot"). Generally, the location of the poles 
varies smoothly with t and continuous loci of roots are readily apparent. Finding _Ra and T{A) follows as indicated in figure 

H 



5 RESULTS FROM THE COMPLEX ANALYSIS 

Root plots were calculated for a range of angles ^ 9 ^ tt. Since the root plots depend upon |A| they are invariant under 
9 ^ 9 — n and this coverage suffices for all possible top-hat models. For the results of the full survey in 9 see Appendix [P] 
The theoretical radius of convergence RA{t) and time of validity T{A) follow directly. 

This section analyses the theoretical convergence for specific open and closed models derived from the root plots. These 
estimates are compared to the time of validity inferred by numerical evaluation of the LPT series. The range of models with 
limited LPT convergence is characterised. The concept of mirror models is introduced to elucidate a number of interconnections 
between open and closed convergence. The physical interpretation of roots introduced by the complexification of the equations 
but lying outside the physical range are discussed. Finally, the special case where the background and the perturbation have 
the same big bang time is analysed. 

5.1 Open models 

Figure IgI shows i?A(i) for 9 = 2.82 and initial scale factor ao — 10""^. All A yield expanding open models for this 9; one choice 
corresponds to the model whose LPT series appeared in figurern(A = 0.01, 9 = 2.82, ao — 10~^). The x-axis is loga and is 
equivalent to a measure of time. The y-axis is log|A|, i.e. the distance from the origin to poles in the complex A plane. In 
principle, future evolution may be limited by real or complex roots. The blue solid line and the red dotted line indicate real 
and complex roots of 77 = 27r respectively. The cyan dashed and pink dot-dashed lines indicate the real and complex roots of 
rj — respectively. Future evolution is constrained by real roots (blue and cyan) in this example. 

The time of validity is the first instance when a singularity appears within the disk of radius A in the complex A plane. 
For the specific case, starting at ordinate A — 10^^, one moves horizontally to the right to intersect the blue line and then 
vertically down to read off the scale factor Oi, — a[T(A)] = 0.179. The time of validity inferred from the root plot agrees 
quantitatively with the numerical results in figure IT] 

AppendixlDlpresents a comprehensive set of results. The time of validity is finite for any open model. As expected, smaller 
amplitudes imply longer times of validity. The poles do not correspond to collapse singularities reached in the course of normal 
physical evolution since the open models do not have any real future singularities. A hint of an explanation is already present, 
however. The green dashed line is 5„ = 1 (or A = 1/sinS) at which point the root switches from rj = 2tt below to above. 
Such a switch might occur if varying the initial velocity transposes an expanding closed model into a contracting closed model. 
But it is expected to occur at (5„ = — 1 not 1. The open models are apparently sensitive to past and future singularities in 
closed models with initial conditions that are transformed in a particular manner. |5.3| explores this interpretation in detail. 
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Figure 6. /?a for 9 = 2.82 and ao = 10"'' (vertical dashed line). To determine the time of validity for LPT expansion with a given A, 
move horizontally to the right of a = ao following the dashed line with arrow and locate the first coloured line with ordinate equal to A 
and then move vertically down to read off the scale factor at the time of validity a^. The specific case illustrated (A = 10~^) matches 
that of the model with problematic convergence in figure IT] The time of validity is correctly predicted. The meaning of the colours is 
discussed in the text. Coloured version of the figure is available online. 
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Figure 7. R/^ for 9 = 0.44 and ap = 10"'^. The line A^; = separates open and closed models. The scale factor at the time of validity 
is a^. For closed models the scale factor at time of collapse is a^. Blue solid line and red small dashed line denote real and complex roots 
of r; = 27r, respectively. The cyan dashed lines denotes the real roots of r; = 0. When the first singularity encountered is real, a^ = a^, 
the time of validity is the future time of collapse. However, when the singularity is complex the time of validity is less than the actual 
collapse time. In the range A^c < A < l\E=Ot there are closed models with a^ < a^. 

5.2 Closed models 

Figure It] presents RA{t) for models with 9 — 0.44 and ao — 10~^. There are several new features. Over the angular range 
9c<0<6t the cosmology is closed for small A (see shaded region in figure M near A = 0). Conversely, a straight line drawn 
from A = within this angular range must eventually cross the parabola E = Q except for the special case = 0. Since the 
velocity contribution to energy E oc A^ while the density contribution oc —A it is clear that eventually S > as A increases. 
The critical value, A_b=o, is a function of 9. Below the brown horizontal dot-dashed line in figure IT] the models are closed, 
above they are open (line labelled A = A_b=o). 

The root plot has, as before, blue solid and red dotted lines denoting the distance to real and complex A poles, respectively, 
for r] = 2ti. The cyan dashed line denotes real roots for r) = Q and does not restrict future evolution. 

For small A real roots determine the time of validity. These roots correspond exactly to the model's collapse time. In 
other words, the time of validity is determined by the future singularity. For example, for A = 0.01, the root plot predicts 
that a series expansion should be valid until the collapse at a = 5.5 denoted by "flu = flc" on the x-axis . This prediction 
is confirmed in the left hand panel of figure |8] The root diagram is consistent with the qualitative expectation that small 
overdensities should have long times of validity because collapse times are long: liruA^-o r(A) — >■ oo. 

As A increases from very small values, i.e. successively larger initial density perturbations, the collapse time decreases. 
Eventually the velocity perturbation becomes important so that at A = A^-c a minimum in the collapse time is reached. For 
A_E=o > A > Arc the collapse time increases while the model remains closed. As A — )■ A_b=o the collapse time becomes 
infinite and the model becomes critical. All models with A > A_b=o are open. 
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Figure 8. The exact solution (black, dashed) and LPT expansions of successively higher order (blue) for two expanding, closed models 
with 6 = 0.44. The left hand panel has A = 0.01. LPT converges to the exact solution at all times up to the singularity at a = 5.5. The 
right hand panel has A = 0.2. LPT does not converge beyond a = 0.38. 



The root diagram shows that for A > Arc, the time of validity is determined by complex not real A for rj = 27r. Closed 
models with Arc < A < A_b=o have a time of validity less than the model collapse time. For example, for A — 0.2, the collapse 
occurs at a = 0.94 but convergence is limited to a ^ 0.38. This prediction is verified in the right panel of figure [8] 

The convergence of LPT expansions for some closed models is limited to times well before the future singularity. This 
general behaviour is observed for 9^ < 9 < 9^ and A,.c < A < A_b=o where both Arc and A_b=o are functions of 9. Appendix 
Id] provides additional details. 



5.3 Mirror models, real and complex roots 

The parametrization of the perturbation in terms of A > and —tt < 9 ^ tt and the complexification of A — >■ A can give rise 
to poles anywhere in the complex A space. When Ra is determined by a pole along the real positive axis, a clear interpretation 
is possible: the future singularity of the real physical model exerts a dominant influence on convergence. LPT expansions for 
closed models with A < Arc are limited by the future collapse of the model and are straightforward to interpret. 

The meaning of real roots for open models is less clear cut. The roots determining _Ra at large t are negative real and 
small in magnitude. Negative A lies outside the parameter range for physical perturbations taken to be A > 0. Nonetheless 
the mapping (A, 6) — >■ (—A, 6* ± n) preserves {5, 5„) and the original equations of motion. The poles of the models with 
parameters (A, 9) and (A, 9 zLn) are negatives of each other. Let us call these "mirror models" of each other. 

For infinitesimal A if the original model is open then the mirror model is closed. Figure [2] shows that the A_b=o line has 
some curvature (in fact, it is a parabola) whereas the mirror mapping is an exact inversion through A = 0. Small A points 
are mapped between open and closed; large A points may connect open models to other open models. 

If the original model is open with limiting pole which is negative real of small magnitude then it corresponds to a future 
singularity of the closed mirror model. For example, the closed model with parameters (A = 0.01, 9 — 0.44) in the left panel 
of figure Is] and the open model with parameters (A = 0.01, 9 = 0.44 — vr) shown in the left panel of figure |9] are mirrors. The 
time of the validity of the open model equals the time to collapse of its closed mirror. 

The notion of mirror models explains other features of the root diagrams. The time of validity of open models was 
previously discussed using figure |6](0 = 2.82). The blue solid line indicated real roots. Such roots are the future singularities 
of closed mirror models lying in the fourth quadrant along 9 = 2.82 — tt = —0.32. As A increases the sequence of mirror 
models crosses the 5v = ~1 line (the horizontal dashed line) to become initially contracting cosmologies and, in our labelling, 
the future singularity switches from rj = 2n to rj = 0. This explains the switch in root label from blue solid to cyan dashed 
seen in figure [6] which occurs at (5„ = 1 in the original model. 

The symmetry of the mirroring is not limited to cases when A is real. It applies for complex A, too. For example, the 
models in the right panels of figures IS] and |9] are mirrors of each other. Their time of validity is the same and determined by 
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Figure 9. Mirror models of the closed models of figure Is] Each graph shows the exact solution (black, dashed) and the LPT expansion 
to successively higher orders (blue) of one mirror model. The original model and the mirror have the same time of validity for the LPT 
expansion. 

complex roots which are negatives of each other. These singularities are non-physical and have no interpretation in terms of 
the collapse of any model yet they limit the LPT convergence in the same way. 

Figure [To] shows the areas of phase space where complex roots determine the time of validity in light red. The area within 
the parabola (light blue) contains closed models. Most of the light blue region has a time of validity determined by real roots, 
i.e. the time to the future singularity. The area with both light blue and red shading encompasses closed models with the 
unexpected feature that the time of validity is less than the time to collapse. 

The area outside the parabola contains open models. The time of validity of the unshaded region is determined by 



real roots. The original observation of LPT's non-convergence for an underdensity (Sahni & Shandarin||1996l is an example 



that falls in this region. For small amplitude perturbations the time of validity is simply related by mirror symmetry to the 
occurrence of future singularities of closed models. The right hand plot in figure |9] is an example of an open model with time 
of validity controlled by complex roots (red shading outside the parabola). 

Finally, some open models (especially those with large A) have mirrors that are open models. Figure ITT] shows mirror 
models {A = 2, 6 = 177r/36) and (A = 2, ^ = 177r/36 — tt). These are initially expanding and contracting solutions respectively. 
The root plot in figure 12 predicts that the series is valid until a„ — 0.0016. The real root with t; = (cyan line) sets the time 



of validity and corresponds to the bang time (the future singularity) of the initially contracting model. 
In all cases, the analysis correctly predicts the convergence of the LPT series. 

5.4 Zeldovich and equal bang time models 

The large expanse of phase space shaded light red in figure [To] suggests that complex roots should play a ubiquitous role 
in LPT applications but the situation is somewhat more subtle. For good physical reasons purely gravitational cosmological 
calculations often start with expanding, small amplitude, growing modes at a finite time after the big bang. The absence 
of decaying modes implies that the linearized perturbations decrease in the past r\ A non-linear version of this condition is 
that the perturbation amplitude is exactly zero at i = 0. The same condition can be formulated as "the background and the 
perturbation have the same big bang time" or "the ages of the perturbation and the background are identical." The condition 



1 
Ho 



y^l 



dy 



=0 [(l + Acose)j/-i + £(A,e)] 



1/2 



2 

3Ho' 



(39) 



^ Our analysis is restricted to the case of initially expanding models, i.e. near A = 0. For initially contracting closed models, similar 
physical arguments motivate a consideration of the behaviour near the initial singularity (not the future bang time). For initially 



contracting open models the epoch of interest is t - 
in the text. 



-oo. These models have large A and are not described by the linear limit discussed 
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Figure 10. The red shaded region denotes part of phase space where complex roots play a role. The solid blue line represents the initial 
conditions which correspond to the background and perturbation having the same big bang time. The black solid parabola separates the 
closed and open models. Coloured version online. 
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Figure 11. Two open models which are mirrors of each other. Each plot shows the exact solution (black, dashed) and the LPT series 
expansion to successively higher orders (blue). The left panel is an initially expanding, open model whose convergence is limited to scale 
factors less than a^ = 0.0016 (arrow). The right panel shows the initially contracting mirror model whose bang time at a„ = 0.0016 is 
responsible for the limitation. 



This is a nonlinear relationship between the two initial parameters A and 9 which is shown by a thick blue line on the phase 
space diagram in figure [lO] We have adopted the name "Zeldovich" initial conditions for the top-hat models that satisfy the 
equal bang time relation. There are a variety of definitions for Zeldovich initial conditions given in the literature. Generally, 
these agree at linear order. This one has the virtue that it is simple and easy to interpret. Note that the blue curve does not 
intersect the region of phase space where complex roots occur except, possibly, near A = 0. 



In the limit of small A eq. ( 39 1 becomes 
A(3sin6'-cos6l) = 0. 



(40) 
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Figure 12. Ha for S = 177r/36 and ao = 10~^. The blue solid and cyan dashed lines denoted real roots with rj = 2n and rj = respectively. 
For A = 2, the time of validity is set by the root with rj = 0, which is the bang time of the mirror model with 6 = 177r/36 — it. See figure 
llll for the evolution of both models. 

The solutions are 9 = 9z± where 6z+ = 2.82 and 6z- = tv — 6z+ = —0.32. The second quadrant solution 6z+ corresponds 
to open models while its mirror in the fourth quadrant 8z- to closed models. Only when A — > can complex roots approach 
the loci of Zeldovich initial conditions but they intersect only in the degenerate limit. 

In the next section, we will show that points starting close to the Zeldovich curve continue to stay near it as they move 
through phase space. Such models have real, not complex, roots. This implies that closed systems along the curve always have 
a convergent series solution. Hitherto, LPT convergence has been studied only for initial conditions close to the Zeldovich 
curve. This is why problems have been noted only in the case of voids. The existence of the complex roots is a new finding. 
All of the above is based on the spherical top-hat model which has a uniform density. 

As emphasized above, there are good physical motivations for adopting Zeldovich-type initial conditions. The fact that 
cosmological initial conditions must also be inhomogeneous (i.e. Gaussian random fluctuations) is not captured by the top- 
hat model. One can imagine two extreme limiting cases for how the simple picture of top-hat evolution is modified. If each 
point in space evolves independently as a spherical perturbation then at any given time one expects to find a distribution of 
points along the Zeldovich curve. As time progresses this distribution moves such that the underdense points cluster around 
the attracting point (—1,0.5) and overdense points move towards collapse. The distribution of initial density and velocity 
perturbations yields a cloud of points in phase space but complex roots never play a role because nothing displaces individual 
points from the Zeldovich curve. Each moves at its own pace but stays near the curve. Alternatively, it is well known that tidal 
forces couple the collapse of nearby points. These interactions amplify the initial inhomogeneities leading to the formation 
of pancakes and filaments. As time progresses motions transverse to the Zeldovich curve will grow. If these deviations are 
sufficient they may push some points into areas with complex roots. In a subsequent paper, we will explore these issues for 
general inhomogenous initial conditions. 



6 LPT RE-EXPANSION 

To overcome the constraints above, an iterative stepping scheme that respects the time of validity is developed for LPT. The 
initial parameters at the first step determine the solution for some finite step size. The output at the end of the first step 
determines the input parameter values for the next step and so on. 



6.1 The Algorithm 

Choose the background (ao. Ho, ^o ~ 1, Yo) and the perturbation {bo = ao, -ffpo, f^po, -'^o) at initial time to- The perturbed 
model is fully characterised by Hpo and Qpo or by So = Ppo/po — 1 and (5„,o ~ Hpo/Ho — 1 or by Ao and 6o. Extra subscripts 
have been added to label steps. 

LPT converges for times t < T{Ao,Oo). Use LPT to move forward to time f, satisfying to < t* < T{Ao,9o). At i*, 
the background and perturbed scale factors and time derivatives are a«, 6*, d*, and b«. The fractional density and velocity 
perturbations with respect to the background are 



S, = (1 + So) 



- 1 



(41) 
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5„,, = -!^-l. (42) 

Re-expand the perturbation around the background model as follows. First, let the time and Lagrangian coordinate for the 
background (inner edge of the unperturbed sphere) be continuous: ii — t, and Yi = Yq. These imply ai — a* and di — a,, 
i.e. the scale factor and Hubble constant for the background are continuous. 

At the beginning of the first step we assumed ao — bo. This is no longer true at the end of the first step. Define a new 
Lagrangian coordinate Xi = Xo6,/a,, new scale factor fei — a^, and new scale factor derivative foi — b^af/b^. These definitions 
leave the physical edge of the sphere and its velocity unaltered 

^physical,* ~ b^Xg = 61X1 (43) 

r physical,, = bt,Xo = biXi. (44) 

The re-definitions relabel the fluid elements with a new set of Lagrangian coordinates and re-scale the scale factor. The 
perturbation parameters are unchanged 5i = S, and 5v,i — (5„,, because physical quantities are unmodified. Consequently, 
Ai = A. and 61 =0*. 

6.2 Flow dynamics in the phase space 

To examine how Lagrangian re-expansion works consider how the Lagrangian parameters A and 6 would vary if they were 
evaluated at successive times over the course of a specific cosmological history. Let 5{t) and 5„ (i) be defined via eq. S and 
apply the second-order equations of motion eqs. ([9| and |10 1 to derive the coupled first-order system 



§ - -f^4l + ^) (45) 

^ = ^{(l + 5„)(l-250-(l + 5)} (46) 

where all occurrences of 5 and S-u are functions of time. From 5{t) and (5„(i) one infers the parameters, A(i) and 9{t). These 
have the following simple interpretation: a Lagrangian treatment starting at time t' has A = A(f') and 6 = 6{t') in the LPT 
series. 

Since the system is autonomous it reduces to a simple flow in phase space. The flow has three fixed points at {S, Sv) = (0, 0), 
the unperturbed, background model, (—1, —1), a vacuum static model, and (—1, 0.5), a vacuum expanding model. Linearizing 
around (0, 0) shows it is a saddle fixed point. The tangent to the _E = curve at the origin is the attracting direction and 
the tangent to the equal big bang curve is the repelling direction. The fixed point at (—1, 0.5) is a degenerate attracting node 
and that at (—1,-1) is an unstable node. The flow vectors are plotted in the left panel of figure [Ts] The blue shaded region 
indicates closed models and red shaded region indicates models where complex roots limit the time of validity for LPT. 

Note that the flow lines smoothly cover the whole phase space. The interpretation is that the continuous relabelling of 
Lagrangian coordinates and re-scaling of the scale factor has the potential to overcome the convergence limitations discussed 
thus far. Otherwise one might have seen ill-defined or incomplete flows or flows that were confined to a given region. 

6.2.1 Asymptotic limits of open and closed m,odels 

The right panel of figure [13] zooms in on the area near the origin. Initial points that correspond to open models starting near 
the origin approach the Zeldovich curve and asymptotically converge to the strong attractor at [5, 5„) = (—1, 0.5). 



Closed models collapse and the density 5 — > oo. In the asymptotic limit, the solution to (46 1 is given hy 5 ^ 5^ -\- K with 
integration constant K. From figure [Ts} the flow lines of closed models that start in the vicinity of the origin trace a parabolic 
path that is parallel and essentially equivalent to the Zeldovich curve. 

The flow shows where re-expansion is needed. Closed model flow lines that start near the origin never pass through the red 
shaded region where complex roots play a role; the time of validity equals the time to collapse and no re-expansion is needed. 
However, closed models that originate in the red region must be re-expanded. The flow suggests that they eventually move 
into the blue region. So even though a closed model may initially have an LPT series with limited convergence, re-expansion 
makes it possible to move into the part of phase space where a single step suffices to reach collapse. 

6.3 Finite steps and feasibility 

This section and the next examine the feasibility of extending a solution from recombination to today. The results will be 
applied to fully inhomogeneous evolution in future paper. 

Let the asymptotic time of validity for an open model be expressed in dimensionless form x = limt->cx) H{t)T{A{t), 9{t)). 
Here, A — >■ y 5/4 and 6 — >■ tan~^ ( — 1/2) = 2.677 and r(A, 6) is determined by the future time to collapse of the closed mirror 
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Figure 13. The left panel shows streamlines of the flow described by eq. ( |46[ |. The colour coding of the plot is same as figure [To| The 
right panel zooms in on the area near the origin which is where all models are located at sufficiently early times. At late times, open 
models move away from the origin towards the attracting fixed point at ((5, 5„) = (—1,0.5). The attraction to the Zeldovich solution is 
shown for a set of initial conditions (yellow, cyan, green and black lines) that begin near but not on the critical trajectory. Closed models 
move out to infinity along the fixed big bang time curve. Coloured version online. 




o 



-2 -1 

Log,o(a(t)) 

Figure 14. Extending the time of validity of LPT. The first step has Aq = 10~^ and 6 = 2.82 and implies scale factor at the time of 
validity a„ = 0.179. Incrementing by half the allowed step gives initial conditions for the second step (Ai,6i) = (0.91,2.68). Note that 
the new time of validity has increased. 



model. The result is x = 2.62 (numerical results in Appendix^}, the time of vahdity is proportional to the characteristic age 
of the background and individual steps grow larger and larger. 

An example shows that the basic effect can be seen even before the asymptotic regime is achieved. Figure [l4] sketches 
the first two steps where the assumed model parameters at the first step are {Ao,6o) = (0.01,2.82). The scale factor at the 
time of validity is a = 0.179. A step with half the allowed increment in time is taken and the system is re-initialised. The 
re-initialisation implies (Ai,6'i) — (0.91,2.68) or {5i,Sv.i) ~ (—0.82,0.4). Afterwards the new time of validity is larger in this 
example. 

The feasibility of the re-expansion scheme can be examined by evaluating the ratio of the time of validity before (T) and 
after (T') a step 

T' 

FigurefTslshows a evaluated along the continuous flow as a function of scale factor for three different starting initial conditions. 
Since a > 3 at all times, starting at initial time ti the time after N steps is roughly t ~ a^f,; > 3^ti. 

Consider, for example, the number of steps needed to extend an open solution from recombination to today. Let tf (ti) 
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a(t) 

Figure 15. The ratio of successive times of validity (a) vs. a{t). The dashed, dot-dashed and dotted Hnes indicate three initial starting 
points (0.5, 0.5), (0, 1), (—0.2, 0.2) respectively. The ratio converges to about 3.6 and the time of validity increases geometrically with A'^. 

be the final (initial) time of interest where tf/ti ~ af/ui ~ 10*'^. Estimating a = 3 implies N ~ logg 10*'^ ~ 10 steps are 
needed. This numerical result for N is an overestimate and one can do better. It is important to recall that it based on an 
arbitrarily high order expansion which achieves an exact solution. If one is limited to calculations of finite Lagrangian order 
and imposes a maximum numerical error at the end of the calculation then more than A'^ steps may be required. At least A'' 
steps are needed for series convergence and more than A'^ steps may be needed for error control. 

One can extend any open model to an arbitrary future time while respecting the time of validity of the LPT series. The 
number of steps is governed by a geometric progression. 

One can also extend any closed model to the future singularity while respecting the time of validity of the LPT series. 



Only a single step is needed for a closed model when the root is real (blue shaded region of figure 13 1. When it is complex 
(the region shaded both blue and red) the model fiows first toward the node at (0, 0) (A decreases) and ultimately reaches the 
region of real roots. Multiple steps will generally be necessary to escape the region of complex roots. An approximate fit (eq. 
(D2|) shows that x ~ T{A, 6)H{t) oc A^ for small A where /3 < —2.5. Both x smd the time of validity increase as the node 



is approached. Time advances at least as quickly as a geometric progression and this is analogous to the manner in which 
the open model steps towards its limit point. However, unlike the open case, once the trajectory crosses into the blue region 
(assuming it does not lie exactly on the unstable attracting trajectory) a single final step is needed. The specific number of 
steps will depend upon the starting initial conditions but will be small because of the property of geometric progression. 

6.4 Demonstrative examples 

LPT re-expansion can solve the problematic convergence in previously analysed open and closed models. 

Open models have asymptotic values of A and 8 and simple evolution. The first section below includes numerical results 
that provide a practical demonstration of the success of LPT re-expansion in this case. Convergence as Lagrangian order 
increases and/or time step size decreases is observed qualitatively. 

Closed models have a somewhat more complex behaviour (before and after turnaround). The second section provides 
both a qualitative and quantitative discussion of convergence. The scaling of the leading order error and the time step control 
which are derived are of general applicability. 



6.4-1 Open model 

Figure 16 investigates the effect of time step and order on the evolution of the open model introduced in figure [I] (A — 0.01, 
6 — 2.82, ao — 10""^). The series convergence breaks down at a = 0.179. The left panel shows an attempt to take a single 
step to a = 1 using successively higher LPT series orders. As expected, higher order terms do not improve the accuracy of 
the description because the time of validity is violated. The middle panel employs three steps to reach a = 1, each respecting 
the time of validity. Now the LPT series with higher order improves the accuracy just as one desires. The right panel employs 
six steps to reach a = 1, each respecting the time of validity. Again, higher order improves the description. Note that more 
frequent re-expansion, i.e. smaller steps in time, improve the errors at fixed LPT order. 
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Figure 16. LPT re-expansion of an open model with Aq = 0.01 and 9o = 2.82. Tlie top tiiree figures show the scale factor for the same 
initial conditions calculated with one step (left), three steps (middle) and five steps (right). The black dots indicate the position of the 
time steps. In the middle and right panels, the solution was advanced 9/10 and 1/2 the allowed time of validity, respectively. The bottom 
figures show the errors for all LPT approximations to b{t) including the unphysical negative ones. The order of the LPT expansion are 
colour-coded according the top left figure. The single step expansion does not respect the time of validity whereas both the three and 
six step examples do. The original expansion does not converge over the full time range whereas the re-expansions do. Coloured version 
online. 



6.4-2 Closed model 



Figure 17 investigates the closed model introduced in figurcl8|(A = 0.2, 6 — 0.44, ao — 10~^). The time of validity is determined 
by a complex root. The first panel shows that the series begins to diverge at a = 0.38 well before the collapse singularity is 
reached at a = 0.94. 

A single time step less than the time of validity is guaranteed to converge as the order of the Lagrangian expansion 
increases. LPT re-expansion utilises a set of such time steps each of which is likewise guaranteed to converge. However, since 
a calculation of infinite order is never achieved in practice, it is worth characterising how convergence depends upon two 
calculational choices one has at hand, the time step and the order of the Lagrangian expansion. 

A single small step beginning at i = io and ending at f/ has leading order error for the m-th order Lagrangian approxi- 
mation PI oc (tf/to — i^^+^^^+i^ where A is the value at the initial time. If the same small interval is covered in A'^ smaller 
steps, the error after N steps scales as N~"^(tf/to — i^^+^^^^+i (^gge Appendix IeI for details). If the step size increases in 
a geometric sequence such that 5t/t is a constant for each intermediate step, then tf — to{l + St/t)^ and the error after N 
steps scales as N{tf/to — l)(5i/i)'"^^A'"^^. This leads to the interpretation that the error per intermediate step scales as 
(5i/t)'"+^A'""''\ Define e = {5t/t)A. The leading order error scales as e™"''^ which is numerically small if e < 1. The sum of 
all the missing higher order terms is finite if St < T, i.e respects the time of validity. 

In a practical application, the initial and final times are not close. A reasonable time step criterion is to choose e < 1 
fixed throughout the evolution and to infer St for a given A. Other choices are possible but St must always be less than the 
time of validity. If e is held fixed throughout the evolution, then the net error after A'^ steps for the m-th order approximation 
ex e'^+^N. 

The number of steps required to go from the initial to the final time can be estimated. As a special case assume that A 



^ Typically, the numerical coefficient is of order unity and varies with m as well as the particular value of 9. For the purposes of a 
discussion of the scaling of the error term, we assume the numerical coefficients to be constant as m and 9 vary. 
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is constant. The time step criterion implies that the number of steps to move from the initial time i = to to the final time 
tf for given e is TV = log(t//to)/log(l + (e/A)). For limited total intervals (f/ — io << to) and small steps (e/A << 1) the 
exact answer reduces to A'' ~ {tf — io)A/e = {tf — to)/(5i. Here St = ef A does not grow appreciably over the interval so the 
estimate for A'^ is a maximum. In this limit, the net error oc e^A. The leading order error for the m-th order Lagrangian 
scheme decreases at least as quickly as e™. 

In more general situations the value of A varies. Once the closed model turns around A increases without bound. For 
fixed e the step size St decreases monotonically to zero as i — > tcoii where tcoii is the time of the future singularity. At any 
order it would take infinitely many steps to follow the solution up until collapse. Consider the problem of tracking the solution 
up to a large, finite value of A = A/. This moment corresponds to a fixed time tf < tcoii in the exact solution. The number 
of steps A'^ < Nmax ~ if/Stf where Stf is the step size for the system near A/; Stf oc e/A/. The leading order error after 
A'^ steps at the m-th Lagrangian order oc e^^^TV < e"^^^ Nmax ~ e^A/. This method of step control forces the leading order 
error at fixed time tf < tcoii to decrease as the Lagrangian order m increases and/or the control parameter e decreases. 

The second and third panels in figure [17] show the runs with e = 0.5 and e = 0.2 respectively. The Lagrangian orders 
are colour-coded; dots show time steps determined by the above criterion. At each order the solution was terminated when 
the numerically determined A > 100 so as to avoid the infinite step regime. This required 32 steps for the e = 0.5 run and 
77 steps for the e = 0.2 run. As the red solid lines illustrates, the first order solution turns around before all other solutions. 
This explains why its step size begins to shrink near the midpoint of the graph. By contrast, all the step sizes for higher order 
solutions are very similar up to that point. 

The numerical errors may be analysed from two points of view. 

(i) A coinparison of different coloured lines (different Lagrangian orders) in a single panel shows that error decreases as m 
increases. This is true in a quantitative as well as qualitative sense. For example, in the second panel at a = 0.64 a plot of 
the log of the absolute error is approximately linear in m, as expected. 

(ii) A comparison of the same coloured lines in the middle and right panels shows that smaller e implies better accuracy. 
Again, this is true in a quantitative as well as qualitative sense. For example, the observed ratio of errors at a = 0.64 for the 
9-th order calculations is 5 x 10~*. To evolve up to this time with e = 0.5 (middle panel) takes 10 steps; with e — 0.2 (right 
panel) it takes 22 steps. The expected ratio of errors is (0.2/0.5)^^^(22/10) ~ 2 x 10"*, the same order of magnitude as the 
observed ratio. 

These comparisons lead to the important conclusion that the leading order error for LPT re-expansion varies with Lagrangian 
order and time step as theoretically expected. 

It is clear that considerable benefit accrues not only from implementing higher order Lagrangian schemes but also by 
limiting time step size (which must always be less than the time of validity). For simple examples like the top-hat it is feasible 
to work to very high Lagrangian order but this is not likely to be true in the context of more complicated, inhomogeneous 
problems. On the other hand, marching forward by many small time steps using LPT re-expansion is generally feasible. In 
the example above the initial perturbation is A = 0.2 whereas a practical calculation starting at recombination would start 
with A ~ 10^"". For the same e the practical application requires more steps for the phase before turnaround but the net 
increase is only a modest logarithmic factor. In fact, most of the steps in the example were taken after turnaround and the 
total number varies with the depth of the collapse. This will continue to be true for the practical calculation. The choice of 
step size and order for such applications will be the subject of a forthcoming paper. 



7 CONCLUSION 

We have investigated the time of validity of Lagrangian perturbation theory for spherical top-hat cosmologies with general 
initial conditions. Using techniques from complex analysis we showed that the time of validity is always limited for open 
models. We also discovered a class of closed models whose time of validity is less than their time to collapse. We introduced 
the concept of the mirror model and derived a symmetry principle for the time of validity of mirror models. For small initial 
perturbations the time of validity of LPT series expansion of an open model corresponds to the collapse time of a closed 
mirror model. 

A qualitative analogy is useful. A single LPT series expansion is similar to a single step in a finite difference approximation 
for advancing a hyperbolic partial differential equation like the wave equation. The time of validity of the LPT expansion is 
analogous to the Courant condition which guarantees stability. In LPT the constraint is an acceleration-related time-scale; in 
the wave equation it is a sound-crossing time-scale. 

We developed the method of LPT re-expansion which overcomes the limitations intrinsic to a single expansion. We 
demonstrated how to iteratively re-expand the solution so as to link convergent series expressions that extend from initial 
to final times. The time of validity of the expansions set the minimum number of re-expansion steps (~ 10) necessary for 
cosmological simulations starting at recombination and proceeding to the present epoch. Finite as opposed to infinite order 
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Figure 17. LPT re-expansion of a closed solution with A = 0.2, 6 = 0.44. The top three figure show the scale factor calculated with 
a single step (left) and multiple steps with e = 0.5 (middle) and e = 0.2 (right) (refer to text for definition of e). The bottom figures 
show the errors for all LPT approximations to b(t) including the unphysical negative ones. The order of the expansion is colour-coded 
as in the top left figure. The single step expansion does not respect the time of validity whereas both the other cases do. The black dots 
indicate the position of the time steps. The original expansion does not converge over the full time range whereas the re-expansions do. 
Coloured version online. 



Lagrangian expansions required extra steps to achieve given error bounds. We characterised how the leading order numerical 
error for a solution generated by LPT re-expansion varied with the choice of Lagrangian order and of time step size. We 
provided a recipe for time step control for LPT re-expansion based on these results. 

Our long-term goal and motivation for this study is to develop a numerical implementation of LPT re-expansion for fully 
inhomogeneous cosmological simulation. Top-hats with Zeldovich initial conditions have special properties with respect to 
LPT convergence. We found that all underdense models must be treated by re-expansion while none of the overdense ones 
need be. However, during the course of an inhomogeneous simulation the density and irrotational velocity perturbations (with 
respect to a homogeneous background cosmology) at an arbitrary point will generally not fall on the top-hat's Zeldovich curve. 
Hence, the convergence of LPT in inhomogeneous applications must be guided by the analysis of more general models. Top- 
hats with arbitrary initial conditions are the simplest possibility and constitute the main focus in this paper. The limitations 
on LPT convergence which we have elucidated in this generic case are considerably more complicated than in the top-hat 
with Zeldovich initial conditions. Our plan is to use the generic time of validity criterion to determine the time-stepping for 
inhomogeneous evolution. This should allow us to develop high-precision simulations with well-defined control of errors. The 
practical impact of a refined treatment of LPT convergence is not yet clear. 

The convergence issues we have dealt with should not be confused with the breakdown when orbit crossing takes place and 
the Jacobian of the transformation from Lagrangian to physical coordinates becomes singular. At that time the flow becomes 
multi-streamed and much of the simplicity and advantage of the Lagrangian approach vanishes. The aim of the current work 
is to make sure it is possible to reach the epoch of multi-streamed flow but offers nothing new on how to proceed beyond it. In 
fact, it may be necessary to include an effective pressure term in the equations to account for the velocity dispersion induced by 
orbit crossing ( Adler fc Buchert|1999 Buchert, Dominguez fc J.Perez-Mercader|1999 1 or to adopt alternative approximations 
for the basic dynamics (such as the adhesion approximation; see Sahni fc Coles|[l995 for a review and references therein) to 
make progress. 
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APPENDIX A: FORMAL SET-UP OF THE SPHERICAL TOP-HAT 

We intend to study an inhomogeneous universe. It contains a single, compensated spherical perturbation evolving in a back- 
ground cosmology. To describe two spatially distinct pieces of the inhomogeneous universe (the background and the central 
perturbation) we invoke the language of homogeneous cosmology. 

Al Description of the background 

The origin of the coordinate system is the centre of the sphere. The background system at the initial time to is set by the 
physical size of the inner edge ri,fi, the velocity fi,_o and density parameter SIq. The Lagrangian coordinate system is extended 
linearly throughout space once the Lagrangian coordinate of the inner edge is fixed. Let the Lagrangian coordinate of the 
inner edge be 

y^*^. (Al) 

ao 
Either choose the initial background scale factor ao and determine the coordinate system or, alternatively, fix Y and infer 
the background scale factor. In either case, the scale factor embodies the gauge freedom associated with the radial coordinate 
system. 

The future evolution of the inner edge of the background is given by rb{t) = a{t)Y. The velocity at the initial time 
satisfies rbfi = aoY. The density at any later time is 

3 

Pb{t) = —^' (A2) 

and the Hubble parameter for the background is 
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Ho = ^ = ^. (A3) 

rb.Q ao 

The evolution of the scale factor is 

a _ AnGpboaf) _ 1 Hgalrio ,.. 

The quantities, r(,_o, '^6,0 , ^o and to along with the choice of the coordinate system, completely specify the background universe. 

A2 Description of the innermost perturbation 

The perturbation can be described by four physical quantities: the physical position Vpfi and velocity Vp^ of the edge (or the 
ratio Hop = rp,o/'"p,o), the density parameter fipo at the initial time to- The Lagrangian coordinate system for the perturbation 
is 

X-^y (A5) 

o(io) 

It can be linearly extended throughout space. 

Like ao, fe(io) embodies the gauge freedom associated with the choice of the coordinate system. Without loss of generality, 
one can pick this gauge to satisfy 

fo(io) = ao. (A6) 

Note that the Lagrangian coordinate systems for the background and perturbation are different. 

Let po and pp,o denote the densities of the background and perturbation respectively. Define the perturbation parameters 

S = ^-1 (A7) 

PbO 

S. = §^-1 (A8) 

-no 
giving 

"-=(IW- ^^'^ 

A3 Inhomogeneous model 



Figure |AT] shows how an overdense and underdense innermost sphere may be embedded with compensation in a homogeneous 
background universe. The assumption that the background cosmology evolves like a homogeneous model, fully described 
in terms of its Hubble constant and density, imposes consistency conditions. At the initial instant the "inner edge" of the 
unperturbed background distribution is at physical distance 7-5,0 from the centre of the sphere. The region with r > rbfi will 
evolve like an unperturbed homogeneous cosmology as long as 

(i) the mass within equals the mass that an unperturbed sphere would contain; 

(ii) matter motions within the perturbed region do not overtake the inner edge of the homogeneous region. 



These conditions which are obvious in the Newtonian context have general relativistic analogues ( Landau fc Lifshitz|l975 K 

Next, consider the innermost perturbed spherical region. At the initial time let Vp^ be the "outer edge" of this region. 
The physical properties and evolution of the innermost region are fully described in terms of its Hubble constant and density 
as long as its outer edge does not overtake matter in surrounding shells. While this is obvious in a Newtonian context there 
exists a relativistic analogue I Tolman|1934 Landau fc Lifshitz|1975 I . 



The inhomogeneous model is incomplete without specification of the transition region between the innermost sphere and 
the background. For the background to evolve in an unperturbed fashion the mass within rtfi must be exactly 47rpo''b,o/3. 
There are many ways to satisfy this requirement. For example, when 5 > a simple choice is to place an empty (vacuum) shell 
for rpfi < r < rbfi so that (ppo/po) = {rb,o/'''p,o)^ = {Y/X)"^. The evolution of each matter-filled region proceeds independently 
as long as the trajectories of the inner and outer edges do not cross. When S < 0, a more complicated transition is required. 



For example, one choice is to nest sphere, empty shell and dense shell (see figure Al I so that the mass within rb,o matches 
that of the unperturbed background. In this case ppo^'pfi ~ f portxi for some / < 1 (the remaining fraction 1 — / is placed 
in the dense shell). Varying the specifics of the compensation region while keeping the properties of the sphere fixed leaves 5 
and 5v, as defined above, invariant. 

For fixed 5 and 5„ the solution b{t) is independent of the details of the transition. Nonetheless, variation in /, Vbfi/rpfi 
and Y/ X all go hand-in-hand. Hence, the extent of time that the sphere's evolution may be treated as independent of the 
matter-filled outer regions also varies. A basic premise of this paper is that it is meaningful to determine the limitations 

© RAS, MNRAS 000,[l]{38] 



26 Sharvari Nadkarni- Ghosh and David F. Chernoff 



Background Dq 

Vacui 




''p.O ~\^~^0^ 




Figure Al. A cartoon showing the physical set-up of the problem. 



arising from the convergence of the LPT series independently of limitations associated with crossing of separate matter-filled 
regions. For a given a S and 5„ this separation can be achieved for specific constructions by choosing the radius and (hence 
velocity) of the inner sphere and the energy of the compensating region appropriately. 



A4 Number of degrees of freedom for the innermost sphere 

If the innermost sphere corresponds to an overdensity then the compensating region can be a vacuum as shown in figure |A1[ 
Having picked the co-ordinate system, having selected equal initial times for the background and perturbation (not equal 
bang times but equal times at which we give the background and perturbation values), and required the correct amount of 
mass, only two degrees of freedom remain: S and 5i,. 

To reiterate, the background and the perturbation can have different big bang times. Setting them equal would imply a 
relationship between S and 5^ and leave a single free parameter. 

If the innermost sphere corresponds to an underdensity then the compensating region is not vacuum but a spherical shell. 
In this case, in addition to 5 and S^, one must specify / or, equivalently, r-p_o- But the solution for b{t) is independent of the 
size of the innermost sphere so, again, only two degrees of freedom remain. 



A5 Preventing shell crossing 

There are two sorts of limitations for the solution of b(t). One is the calculation-dependent limitation arising from the 
convergence properties of the Lagrangian series expansion. It involves the scale factors only. The other is a physical limitation 
arising from collisions of the innermost region with surrounding non- vacuum regions (either the background or a compensating 
shell). We show that it is possible to delay the epoch of collisions indefinitely without altering the evolution of the innermost 
region. 

Fix Ho, Hpfi, pa and pp,o- This implies that the expansion parameters in LPT, 5 and 5^, and the time of validity of the 
LPT solution are all fixed. Consider the case of an overdensity surrounded by vacuum. To stave off the collision of the outer 
edge of the innermost region with inner edge of the homogeneous background hold rj,_o fixed and reduce rp,o- The velocity 
■f-pfi = HopTpfi becomes arbitrarily small. The time for the edge to reach any fixed physical distance increases without bound. 
Shell crossings may be put off indefinitely. However, we have altered the mass within the innermost edge of the background 
so we add back a thin, dense shell just inside rbfi and set it on a critical trajectory outward. This accomplishes our goal. 
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The case of the underdensity surrounded by a compensating shell is identical. First, we must make sure that the com- 
pensating shell does not overrun the homogeneous model. Choose the shell to be thin, fix its initial physical distance from 
the centre and adjust is velocity (based on how the interior mass changes) to give a critical solution. The two power laws, 
one for the compensating shell and one for the innermost boundary of the homogeneous model, cannot cross in the future. 
Second, as above, note that reducing rp.o reduces the outward velocity of the edge so that it takes more time to reach the 
initial position of the compensating shell. The time can be made arbitrarily long. 

The limitations in LPT convergence are completely distinct from those associated with physical collisions in inhomoge- 
neous model. 



APPENDIX B: SERIES EXPANSIONS FOR A FUNCTION OF TWO VARIABLES 

In this section we elucidate by example some qualitative features of the expansion of 6(i, A), the central quantity in the 
Lagrangian treatment of the top-hat. We assume a very simple form denoted f{t, A) and look at convergence with respect to 
expansions in t and A. Let 

/(i,A)=i^''^(i+A)'^'. (Bl) 

The series expansion of this function around A = at fixed t is 

/ ^ il/3 ^ ^4/3^ _ 1 7/3^2 AilO/3^3 _ j^^l^/S^^ + ^^16/3^5 ^,^6) (g2) 

•' 9 81 243 729 ^ -^ ^ ' 

which is supposed to mimic the Lagrangian expansion in A. One can also expand the function as a series in t around t = ti 



1 ,/3 ^ (2Ai, + l)(i - U) ^ (-A^i? - AU - 1) (t - Uf 



S^A+j^yhl^^'iAU + l)^ 

Both expansions involve the complex power z^'^. There are two branch cuts which extend to 2; = so at A = —1/^ the 
function is not analytic. Additionally, the expansion in t is not analytic at f = 0. 

The efficacy of various expansions are illustrated in figure |B1[ In all the plots the black dotted line indicates the exact 
function. The top left panel shows successively higher order series approximations in A as a function of t for the specific case 
A = 1/10. The question here is whether the pole at a given time lies with a disk of radius 1/10? The location of the pole is 
A = — 1/i so the answer is "yes" when t > 10. This pole interferes with the convergence of the series expansion for A — 1/10. 
The figure demonstrates the (future) time of validity is i < 10. 

The top right panel shows the series in A at a fixed t = 1/10. The question here is how big a perturbation will converge 
at i = 1/10? Since the location of the pole is A = — 1/i the radius of convergence at the indicated time is 10. Perturbations 
with |A| > 10 are not expected to converge and the figure shows that this is indeed the case. 

The bottom left panel shows the series in t expanded around ti — 2 for fixed A — 1/10. The poles are at f = —10 and 
t = in the complex t plane. The expected radius of convergence is min(|2 — 0|, |2 — (— 10)|) = 2 or ii — 2 < i < ti -(- 2. As 
seen in the plot, the series converges only in the expected range (0, 4) 

The bottom right panel shows the series in t expanded around ti — 2 for A = —1/3. The poles are at f = 3 and i = in 
the complex t plane. The expected radius of convergence is min(|2 — 0|, |2 — 3|) = 1 or ti — 1 < i < ii + 1. As seen in the plot, 
the series converges only in the expected range (1,3). 



APPENDIX C: PARAMETRIC SOLUTION 

The background model has scale factor ao and Hubble constant Hq — ao/ao- The model, perturbed in density and velocity, 
is parameterized by A and and has scale factor b{t) . For the choice of coordinate system given in the text the second order 
equation for b is 

b lHial(l + Acos6) 



b 2 63 



(CI) 



with the initial conditions that at t = to, 6(to) = ao, 6(fo) = fflo(l + Asin^). The scale factor ao and the velocity of the 
background do at the initial time to are positive. The parametrization of 6(fo) allows either positive or negative values where 
A is non-negative and — tt < ^ n. The quantity (1 + Acos6'), proportional to total density, is non negative. 
This equation once integrated is 
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Figure Bl. Series expansions in t and A for an illustrative function /(t, A) (see text). The black dotted line indicates the exact function 
/ and the blue solid lines indicate successive approximations. The top left and right panels are series expansions in A around A = 
plotted as a function of t (for A = 1/10) and function of A (for t = 1/10) respectively. The bottom left and right panels are series 
expansions in the t around t = 2 plotted as functions of t for A = —1/10 and A = —1/3 respectively. 






(l + Acos6i) (1 + Asin^ 



(l + AcosS) 



ao 



The combination 



E{A,e) = (1 + Asin6')^- (1 + Acose*) 



(C2) 



(C3) 



is proportional to the total energy and determines the fate of the system. If E{A, 6) > 0, the model is open and if ^(A, 9) < 0, 
the model is closed and will re-collapse eventually. Four cases (positive and negative E, positive and negative 6o) are shown 
in figure [5] 



CI Initially Expanding Solutions 

The expanding case with 60 > for open models [E > 0) has solution 
ao (1 -I- Acos^ 



2 E{A,9) 
1 (1 -I- A cos 



(coshr; — 1) 
^(sinh77-r,)+t+„3(A,e) 



biv,A,9) 

*(^'^'^) - 2Ho E{A,9)V^ 
and the singularity 6 = occurs at r; = 0. For closed models [E < 0) the solution is 
ao (1 + Acosf 



h{'n,A,9) 



2 \E{A,9)\ 



-(1 — cosry) 



(C4) 
(C5) 

(C6) 
(C7) 



/ » ^^ 1 (1 + A cos 6), . , + , ^ „-, 

tiv,A,9) = 2/^),^(A;^('/--.)+t+_(A,e). 

For closed models, the convention adopted sets 77 = at the singularity nearest in time to to- For both models, the time at 
7; = is denoted tf^„„- For closed models the time at 77 = 27r is denoted t^^n- 

At the initial time the solutions (both open and closed) satisfy 6(io) ~ 0.0, b{to) — do(l + Asin6l) and t = io- The condition 
b{to) — ao sets the value of the parameter at the initial time r;o. The velocity condition is then manifestly satisfied from the 



form of eq. ( C2 1. The condition t = to at 77 = 770 sets the value of the bang time 
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bang 



to- 



^-^ -[smhrio ~rjo) E > 0. 



{ 2Ho £;(A, 9)3/2 

The bang time for the model can also be written as 
"" db 



bang 



= to 



t=o (&2)(i/2'' 



(C8) 



(C9) 



where b is given by eq. (C2| with the sign for the square root positive. The age of the model since its birth is 

u,4A,9)= f "-^-^^^z ;:r:' ,:: ■ (cio) 



db _ /■"='"' db/dy ■ dr] 

Uo (i2)(V2) - 7^^„ (fe2(^))(l/2) 



Inserting the appropriate parametric solution, one can verify that the bang times obtained from ( C8 1 and ( C9 1 are identical. 
Generally t+ „^ / 0. 

The velocity at the initial time is 



A - A 1 P|l/2 J l-cos7,o -C' ^ U 

bo - ao\E\ <; ,inhC E>0. 



(Cll) 



cosh T^o ~ -^ 

First, 6o > implies 770 > 0. Second, if the age of the model increases, 77 increases. For the open solution if 77 varies from to 
00 time increases from t^j,„ to 00. For a single cycle of the closed solutions, ?7 increases from to 27r and time increases from 
tt tof+„. 

bang coll 



In summary, the parametric solutions solve eq. (Cll and eq. (C2| for the specified initial conditions. As a final useful 
step, rewrite eq. ( |C9| | by defining y = b/ao 



bang 



to- 



1 



dy 



Hojy^o l{l + Acos9)y-^ + E{A,e)] 



1/2 



(C12) 



which follows from eq. (C2l and uses the same positive square root convention. 



C2 Initially Contracting Solutions 

Next, consider the case 60 < 0. The parametric solution for _B > is 
ao (1 + Acos^ 



b{v,A,9) 

t{n,A,e) 

and for iJ < is 
6(77, A, e) = 

t{n,A,e) = 



2 E{A,e) 

1 (l + Acos 



2Ho S(A, 61)3/2 



(cosh 77 — 1) 

'\-sinh77 + 77) + i- (A,e) 



(C13) 
(C14) 



ao (1 + Acos^ 



-(1 — cos 77) 
^(-77 + sin77) + i-„g(A,e). 



(C15) 
(C16) 



2 \E{A, , 
1 (l + Acos( 
2Ho \E{A,e)\3/2 

Again, for closed models, the convention adopted is that the singularity nearest to to corresponds to 7; = 0. The time at 77 = 
is t^^^ and the collapse time for closed models is f"^;;. 



The parametric form of the solutions satisfies eq. (Cll and eq. (C2|. Just as in the previous case, the initial conditions 



set 770 and t^^^ . Since the singularity at 77 = lies to the future of io, 
db 



bang 



to + 



9 
b=an 



b=o ib^Y'/^^ 



(C17) 



where, 6 is given by eq. (C2|. The sign of the square root is chosen to be positive and the integral is a positive quantity 



which is added to to. For closed models the singularity at 77 = 27r lies to the past of io at f^o„. In this case (see figure l5| the 
labelling implies tj^jj < to < tf^anq- Although this might seem backwards, it facilitates combining the open and closed models 
into one complex function as was done in the positive bo case. The initial velocity is 

bo = ao\E\^^^ • 



1 — cosh 770 
3in'7o 



E>0 
E <0 



^ , ,. (C18) 

cos r]Q~l 

The initial velocity 60 < implies 770 > 0. For the age of the model to increase, 77 must decrease. Conversely, if 77 increases, 
the time in the open model decreases from to to —00 and the time in the closed model decreases from to to t~^ii- 
A table summarising the properties of the physical solutions with t/ = \ri\( — rjC, follows. 
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Closed 


Open 


If rj increases 


If t increases 


ibang f-O 


bo >0 


C = i 


<; = i 


t increases from t'^^„ 


r; increases to oo or 2it 


<0 


bo <0 


C = i 


( = i 


t decreases from fj^^^ 


7] decreases to 


>0 



C3 Analytic Extension of the exact solution in parametric form 



The differential eq. (321 was solved numerically over the range 0i^9^-K,0<S< 100 and — tt < 4> ^ tv where A — Ae"*. For 
each value of (A, (j), 6), the numerical solution matched one of the two possible parametric forms. 
Omitting the explicit functional dependence on A and the following abbreviations are useful 

j = (1 + A cos 6*) 

(1 + Asin6l)^ 
h — - — 

J 

E = (h-l)j. 

The two possible parametric forms that agree with the numerical solution are 
ao j 



Hv) = 



t(r?) = to ± 



where 



2 [-E]^ 
/ 1 



cos ri) 
j 



\2Ho [-E] 



3/2 



(^- 



1?;) - ta 



TJase — 



Ho 



]Vh- 



-r- sinh 

[Ef/^ V J 



(C19) 
(C20) 
(C21) 

(C22) 
(C23) 

(C24) 



The branch cut lies along the negative real axis for all fractional powers and from — ioo to — i and +i to ioo for the inverse 
sinh function. 

The prescription for the correct form is for the choice of the ± sign in t eq. ( |C23[ ) and denoted t+ and t_. The correct 
form depends upon 6, <j), arg[h] (the arg is defined to be between — vr and vr) and the (real) value J = j when (^ = or tt. The 
figure [CT] shows the upper half plane for the perturbation partitioned into areas where the complex extension of the solution 
has one of two forms. The lower half plane has the same structure inverted through the origin. The horizontal red dashed line 
denotes Asin6' = 1 and the vertical red dashed lines denote AcosS = ±1. In some areas a single form applies as marked but 
in the central area both occur. The detailed prescription is 



s; e s; 7r/4 



= 7r,|A| sine < 1 and j < t- 
otherwise ta 



t = < 



7r/4 < 6* < TT < 



< <^< TT 


and arg h > 


t 


— TT < (f) <Q and arg h < 


t 




f cos e > 




4> ~ Q and 


^_ cos 61 < and j > 
( cos 6^ < 


t 


(j) = TV and 


or 

(^ cos > and j < 


t 


otherwise 




t 



(C25) 
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Figure CI. This fig ure de scribes one aspect of the analytic extension of the exact solution. For a given real A, the complex extension 
A — > Ae**^ obeys eq. l|C25l with two possible forms t+ and t_. The choice depends on 4>, A, 9. For some (A, 9) a single form is sufficient 
for all <p; for other values both forms are needed. This figure illustrates how the upper half plane is partitioned based on this property. 
Coloured version online. 



APPENDIX D: NUMERICAL SOLUTIONS 



Dl Algorithm 

The initial conditions are parameterized by A > and — vr < 9 ^ n. The transformation 6 — > tt ± 5 and A — >■ —A leaves the 
solution unchanged. At any time the roots for 9 and -k ±9 are negatives of each other. The root plot only depends upon the 
absolute value of the root so the plots for 9 and -n ±9 are identical. It is sufficient to consider the upper half plane. 

For a given 9 the algorithm to map out R/\{t) is the following: Vary A from to an arbitrarily large value (~ 100) in 
small increments. For each A select A = Ae"* by varying the phase angles 4> over the range to 27r. For each A evaluate 
t{ri) a.t T) — Q and rj — 2tt calculated according to eq. ( |C25| |. Finally, hunt for solutions that set the imaginary part of t to 0. 
This last step involves one-dimensional root-finding in (j) at fixed A. A solution leads to a specific pair (i, A) that is a pole in 
the function b(A,t). 

Roots with t > to limit future evolution; those with t < to limit backwards evolution. Both sets are shown in the results. 
Roots are classified based on whether they are real or complex. For closed models the real roots can represent a singularity 
that is nearby {rj — 0) or far away {rj = 2n) from to- This classification at the initial time is independent of whether the 
singularity is in the past or future and is independent of whether the model is expanding or contracting. For open models the 
real roots are always considered nearby (t; = 0). 

In what follows the numerical answers are first described in qualitative terms. In the next section simple analytic estimates 
for the time of validity are developed. 

Figure [Dl] shows the root plots on a log-log scale. Sixteen panels, each with a particular value of 9 listed at the top, are 
displayed. The x-axis is logj^Q Hot and the y-axis is the log of the distance of the singularity from the origin in the complex A 
plane. The initial time, Hoto = 2/3, is marked by the vertical black dashed line. 

For each 9, the shaded region indicates the range of A that gives rise to closed models. Figure [2] shows that closed models 
occur only for 9 < 9^ = 0.463 in the upper half plane so only some of the root plots have shading and then only at smaller A. 

The colour coding of the dots indicates four types of roots: real and complex roots where rj = 2-n are in blue and red, 
respectively; real and complex roots with t; = are in cyan and pink, respectively. The radius of convergence at the initial 
time io is infinite, i.e. the Lagrangian series is exact at the initial time by construction. At times very close to the initial 
time the root loci lie ofi^ the plot. Only the roots to the right of Hoto are relevant for forward evolution and, conversely, only 
those to the left are relevant for backwards evolution. The discussion is focused on the case of forward evolution but it is 
straightforward to consider the restrictions on marching backwards in time. 



The phase of the root (of smallest magnitude) appears in figure D3 When closed models have real roots they are positive; 
when open models have real roots they are negative. However, some open and closed models also possess complex roots. The 



set of models with complex roots (of smallest magnitude) is evident from the shading in figure 10 The phase of each root of 
smallest magnitude in figure [Dl] is indicated by the colour shading in figure [D3| 

There are horizontal dashed lines with colours green, blue and purple in figures Dl and D3 indicating \5v\ = 1, \5\ — 1 
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Figure Dl. Root plots for 6 in the range ^ S ^ tt. In each plot the abscissa is logjQ Hot and the ordinate is the logarithm of the 
magnitude of the root. The vertical black dashed line marks the initial time. The shaded area corresponds to closed models. The blue 
and red points show real and complex roots with r; = 2n, respectively. The cyan and pink show real and complex roots with r; = 0, 
respectively. The green and purple dashed lines are |<5„| = 1 (A = | sinS|~^) and |5| = 1 (A = | cosS|~'^), respectively. The blue dashed 
line indicates the switch between two forms and a single form of the parametric solution at A = |2secS — cscS|. Coloured version online. 



and the transition between one and two complex forms, respectively. For each 9 the Hnes mark the imphed, special value of 
A. These dashed lines also appear with the same colour coding in figure [02] 



The roots in figure Dl will be analysed in the range < 9 ^ i"/4, 7r/4 < 6 ^ 'r/2, -k/2 < 9 ^ n. 



Dl.l Q^9(^Tx/A 

The top left panel in figure [Dl] has 9 — Q\ the blue dots indicate real roots with r) — 2tt; the blue shading indicates a closed 



model; the phase is positive (top left panel in figure D3l. Only a single branch is evident. In sum, each root is the collapse 



time of a closed, pure density perturbation. For an expanding model, rj = 2tt implies that the root is the future singularity. 
That S = is a special case can be seen by consulting figure [02] a ray starting at A = with 9 = never intersects any of 
the other lines of the diagram. In general, each time a ray crosses one of the lines there is qualitative change in the properties 
of the roots. 
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Figure D2. Several conditions determine the nature of the roots in phase space. The most significant are schematicaUy illustrated here. 
The green horizontal lines are |(5i,| = A| sinS| = 1; purple vertical lines are |5| = A| cos9| = 1; the black curved line is the B = critical 
solution. The red lines Arc mark where real roots associated with closed models (or closed mirror models) transform to complex roots. 
The blue dashed lines mark the division between one and two complex forms (see also figure [Clj l . Physical models lie to the right of 
& = —1. Expanding models lie above 5^ = —1. The intersection S = Sy = 1 occurs at 9 = 7r/4. The point P near 9 = 0.84 is the meeting 
of (Si, = 1 and A^c. Coloured version online. 



For < 6 ^ 7r/4 a great deal more complexity is evident in figure Dl First, consider a ray emanating with small angle 
(tan( 



Q <0 <9J in figure 



D2 



— 1/2 is the slope of the E = Q line at the origin). Eventually such a line will cross the black 
line which is the E = Q critical solution labelled A_e=o- For small A the models are closed; for larger A they are open. In 
figure [dT] this distinction corresponds to the the blue shading (closed models) at small A versus the unshaded (open) models 
at large A. 

Within the shaded region note that two branches of real roots are present beyond a given time; at large t (asymptotically) 
the lower branch is A — )■ and the upper branch is A — )■ A_b=o- The lower branch sets the time of validity for small A. Each 
root is the collapse time of a closed model which has both density and velocity perturbations at the initial time. 

As A increases the time of validity inferred from the lower branch decreases. At the critical point A = A^c, the two real 
branches merge and connect to a branch of complex roots (intersection of red and blue points). For A > Arc, the complex 
roots determine the time of validity even though the upper branch provides a real root. The complex roots do not have a 



© RAS, MNRAS 000,[T]{38] 



34 Sharvari Nadkarni- Ghosh and David F. Chernoff 

direct physical interpretation in terms of future singularities of physical models. On figure [D2] the ray emanating from the 
origin at shallow angle crosses the red dashed line labelled Arc at this critical point. 

Physically, when A exceeds Arc, the velocity perturbation dominates the density perturbation in the sense that the 
collapse time begins to increase. The real root corresponds to the future singularity of the model. As A increases further, the 
solution eventually becomes critical (infinite collapse time). The particular value where this occurs is A_e=o and it corresponds 
on figure [D2l to the ray crossing the labelled black line. Within the entire range Arc < A < A_e=o the complex root determines 
the time of validity. So, even though any model in this range is closed and possesses a real future singularity, the time of 
validity is determined by the complex root. This gives the sliver on figu re|10| which is the overlap of light red and blue shadings. 

Both Arc and A_b=o decrease as 6 ^>- 9f as is evident from figure |D2| and both vanish at 6'^ . On figure Dl the real roots 
completely disappear and only the complex roots are present, i.e. the two real branches have been pushed out to infinite times. 
The panel with 9 = 57r/36 — 0.436 is numerically closest to the critical case 6^ = 0.464 and the real branches are just barely 
visible at the right hand edge. 

For the rest of the upper half plane 9^ < 9 ^ n the ray no longer intersects any closed models. 



For 9^ < 9 < 7v/4 the real roots reappear and move back to the left in figure Dl (see panel with 9 = 7r/6). Now, however. 



the roots are negative (see figure D3l. This is a manifestation of mirror symmetry which relates the negative real roots of 
an open model to the positive real roots of a closed model. At large t the two branches have A — >■ and A — > A^^o and 
are completely analogous to the real branches just discussed for closed models. The separation between the two real branches 
increases as 9 increases and the solution loci shifts upwards in A. And just as before the two branches join and meet a complex 
branch. The second red dashed line Arc in figure [02] shows the real to complex transition for the roots for the open models. 



This behaviour might be expect to continue for 7r/4 < 9 < tv but there is an additional complication: the analytic 



extension involves two forms. As the ray sweeps counterclockwise in figure CI it crosses (5„ = 1 (horizontal dashed line and 
the curved blue line. These are also schematically illustrated in figure [02] 

D1.2 -k/4<9 < -k/2 

All physical models are in this range are open. Real roots have a straightforward interpretation in terms of the mirror models. 
Although some of the analysis described for 9 < 7r/4 continues to apply several additional complications ensue. To understand 
them it is useful refer to the phase space picture shown in figure |D2[ As 9 increases, the point where Arc meets 5„ = 1 is 
labelled P. 

For a fixed 9 consider increasing A from small values near the origin to oo. The order in which this ray intersects the 
green (5„ — 1), purple {5 — 1), red (Arc) and blue (one or two complex forms) curves will correlate with the change in roots. 

The roots are negative real for small A. They correspond to the collapse time of a closed mirror model. Increase A 
and ignore Arc. When the (5„ = 1 line is crossed, the sign of the closed mirror model's velocity switches from expanding to 
contracting. This just means that the labelling of the future singularity switches from further away {rj = 27r) to nearer (r; = 0). 
Now recall A < Arc implies real roots and, by definition, 5„ — Asin^. Hence, Arc(^) > l/sin6l implies that the label switch 
occurs just as outlined. On figure |D2| rays counterclockwise of point P belong to this case. This is responsible for the switch 



from blue (real rj = 27r) to cyan (real ?; = 0) roots at the green line in figure Dl for 9 = 7r/3 and 177r/36. 

Conversely, if Arc(^) < 1/ sm9 the roots are already complex and the label switch occurs between the corresponding 
complex roots. There are no pictured examples in figure iDl] 

In the previous section with the < 9 ^ 7r/4, the physical interpretation of Arc (as A increases) was that the velocity 
contribution to the perturbation became dominant in the original model if the model was closed or in the mirror closed model 
if the original model was open. In latter case the mirror models were initially expanding. Now, the same idea continues to 
apply in the regime 7r/4 < 9 < 0.84. Here the transition from real to complex roots occurs before the Ji, = 1 line is crossed. 
The significance of Arc is that it marks the increasing importance of velocity perturbations in the closed expanding mirror 
models. 

However, for 0.84 < 9 < 7r/2 as A increases the open model crosses Sv = 1, the mirror model swaps from 77 = 27r to 
and the roots (real) corresponds to the real future singularity of a closed, contracting model. As A increases further, first the 
mirror model becomes critical and then an open model contracting to a future singularity. While the magnitude of A grows 
larger than a critical value the velocity perturbation dominates the mirror model dynamics. When A > Arc the roots switch 
from real to complex. At this point the contracting mirror model can be open, closed or critical. 

Note in figure |D2] that Arc asymptotes to the vertical purple line 5 = 1. The corresponding mirror model hits the line 
5 = —1 in the third quadrant. This is the limiting vacuum solution. Although there are no physical models beyond the analytic 
extension continues and the roots change from real to complex. All open models with 9 < 7r/2 see a transition to complex 
roots as the mirror approaches the vacuum solution. 

Finally figure |D2| shows as a blue curve the point at which there is a switch in complex form of the analytic extension. 
Here, the complex roots switch from r] — to rj = 2n. The roots remain complex and since there is no physical interpretation 
and it is irrelevant whether they belong to 77 = or ry = 27r. 
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In figure Dl the panel with = 7r/3 and 177r/36 show these transitions: the blue to cyan transition at the green dashed 
line is the mirror model switch from expanding to contracting; the cyan to pink transition at the purple dashed line is the 
mirror model moving through 5 = —1; the pink to red transition is the switch from two to one complex roots and r; = to 
rj = 2-n. 

D1.3 6 = 1^/2 

At 9 = T^ 12, only real roots of r; = are present for large A. This is a special case in that a ray only intersects one special 
line (5i, = 1 in the upper half plane. 

Dl.i 7r/2 < 6* s; TT 

All models in this range also correspond to open models. Like the previous cases, small A have real, negative roots with 
r\ — 2-K. The mirror models in this case lie in the fourth quadrant. The crossover of real roots from r] = Q to r] = 2ii occurs at 
5v = ^, however, unlike in the earlier case, the line 5 = — 1 is never approached by the mirror models in the fourth quadrant. 
As a result, there is no switch from real to complex roots and all models have real negative roots. The rj — 2n roots for small 
A are collapse times of initially expanding closed mirror models and the r; = are future singularities of initially contracting 
closed and open mirror models for intermediate and large values of A respectively. 

D2 Numerical Results 

Here we present numerical formulas that give the time of validity for any initial A and 6. Real roots occur for small A when 
< < n/2; and they occur for all A when 7r/2 ^ ^ ^ tt or ^ = 0. Real roots correspond to past or future singularities of 
physical models and are known exactly. 



Figure Dl shows that complex roots occur < 9 < n/2. In the range n/A < 9 < ■k/2 figure D3 shows that the phase of 
the complex roots is very close to vr. We can approximate these roots as real, negative roots. Conversely, figure [DS] also shows 
that in the range < 9 ^ 7r/4 the phase is not close to or n. These roots are complex only when A > Arc. First, we fit Arc 
by 

Arc,app{9) = |0.41csc^6'(cos6l-2sine) + 3.57(cos6'-2sin6')(sin6i)'*-^^| . (Dl) 

We cannot approximate the time of validity with the results for physical cases but it turns out that the numerically derived 
time of validity is insensitive to 9 in the range < 9 < n/4 and may be fit 



2 

Hotappi^) ~ ^ + < 



A2B 
1.514 

/\2.82 
1.778 

63095 
2x10*^ 



< A^ 1 

1 < A^ 2 

2 < A ^ 5 (D2) 
5 < A sC 10 

A> 10. 



Using these quantities, the table below gives an approximation to the time of validity, Tapp, for all values of 9 and A. The 
times for collapse and the bang times are equivalent to eq. ( |C23[ | and reproduced here for convenience: 

t.MA,9) = to + ^^l^±^^i2n)-U,4A,9) (D3) 

t,-„,(A,e) = to+ta,4A,9) (D4) 

where, 

,^n^ 1 /t; X ^ /(1 + Asin6')2 1 (1 + A cose) . , _i I E{A,9) ,^^, 

E{A,<i),9) = (l + Asin6i)^-(l + Acos6') (D6) 

The error in the fit is estimated as 

£ = ^^p^. (D7) 

If f > then the approximation is conservative in this sense: the approximate time of validity is less than the true value. 
Conversely, if f < 0, then the approximation overestimates the time of validity. Using the above fits the worst case is £ ~ —0.02. 
We always use a time step St which satisfies St < 0.98Tapp so that the inaccuracy in the approximation is irrelevant. 

© RAS, MNRAS 000,[l]{38] 



36 Sharvari Nadkarni- Ghosh and David F. Chernoff 

Table Dl. Approximation to time of validity, Tapp{A, 9), for ^ 6 ^ tt. Note that A^capp is an approximation to A^c in eq. (JDljI. 



Parameter range 



E{A,e)<o t,M(A,e) 

0<e<7r/4 0<A< A,.e,app E{A, e)>0 tcoiii-A, e) 

o<^<R^ tcou{-A,e) 

^>\ 'TnL::i' \ sR[t,„„(-A,e)] 



0<A=£^j5^ t,M{-A,e) 

7r/2 ^ e s£ TT A > ^i^ fr (-A,e) 

/ ^ ^ I sin 0\ bang ^ ' -* 

APPENDIX E: ERROR CHARACTERISATION OF THE LAGRANGIAN SERIES 

We want to characterise the errors associated with calculating the solution at time tj given some fixed initial conditions at 
time io- Errors arise because any real calculation involves truncating the Lagrangian expansion. We want to compare the errors 
that result from different choices of truncation order and of the number of re-expansion steps assuming all series expansions 
are convergent (i.e. all respect the time of validity). Let Nm represent the final physical coordinate generated with a m-th 
order calculation using A'^ steps. Ultimately, we seek to characterise differences like Nm — N'^,' ■ The quantity loo is the exact 



EO.l Single step 

The Lagrangian series solution for a single step has the form 



loo = a(f)ll + ^-^AnXo, (El) 

where each fe''^ satisfies 

6« - ^^°^ = 5« (E2) 

and initial conditions are specified at i = io- The initial conditions at each order and the forms for the first few S^^' are given 
in the text. 

li tf — to = St « to, then the solutions can be expanded in the small parameter St /to- The solutions are 

b^'\t)/a{t) ~ c(^)J^, (E3) 

to 

b^'\t)/a{t) ~ c<'' f-V fori>2. (E4) 

The coefficients c*-'' depend on the angle 6 and have a weak dependence on the Lagrangian order. The difference between the 
exact answer and the m-th order approximation for a single step is 

\i=m + l / 

As long as tf is within the time of validity of LPT, by definition, the LPT series converges and from the equation above, the 
leading order error scales as ~ (St/to)"^^^ A"^'^^ (order terms first by powers of A and then by powers of St/to). 

E0.2 Multiple steps 

In general for a practical application one is limited to working at a finite Lagrangian order. In such cases, it becomes necessary 
to ask if convergence can be achieved by working at a finite Lagrangian order with increasing number of steps. 

First, we outline the calculation. The initial data is subscripted by "0" . For example, let the initial perturbed scale factor 
be 6o = b{to), the initial background scale factor ao, the initial density contrast So and the initial velocity perturbation SvO- 
The Lagrangian expansion parameter Ao and angle 9o follow from the relations So — Aocos^o and Svo ~ Ao sin 6^0. The 
physical coordinate is ro — boXo; for given ro the initial Lagrangian coordinate Xo is fixed by choosing bo to be equal to ao. 
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Consider taking A^ steps from initial to final time with an 77i-th order Lagrangian expansion. Assume that the final time 
is within the time of validity of the Lagrangian expansion. For definiteness, let the j'-th time be tj = io/3"' where /3 — tf /to 
(so tjv is just the final time tf). This geometric sequence of increasing steps is well-suited for an expanding background with 
a small growing perturbation. The scaling of differences like Nm — ^aa, {N + 1)™ — A'^^ and Nm+i — N,n with TV and m are 
all of interest. We expect the same scaling of these differences with TV and m for any uniformly refined set of time steps. 

A finite order Lagrangian expansion accurate to order m is a truncated representation of the full Lagrangian solution 

m 

6(i) = ^6W(t)A^. (E6) 

At the beginning of the first step the scale factor at to is advanced to ii and written as b{to — >■ ti; Ao, So). Note the explicit 
dependence on the perturbation parameters at io- Abbreviate the scale factor and its derivative for the truncated expression 
as 6 and b. The background scale factor at time ti is ai . At the end of the first step the Lagrangian coordinate Xi and the 
new foi are inferred as described in the main body of the text by re-scaling quantities calculated at ti. The new 61 is not b. 
The net result for the full step to — > fi is 

Xi = Xo— (E7) 

ffli 

&i = ai (E8) 

bi = ^ai (E9) 



<5i 



(l + 5o)(^)'-l (EfO) 

5„i = ^-1. (Efl) 

aib 

The newly defined quantities subscripted by "1" will be used to initiate the next step. The updated perturbations imply new 
Lagrangian expansion parameter and angle according to 

Ai = VW+Wi (E12) 

cos6li = -^ (Ef3) 

Ai 

sin 6)1 = ^. (Ef4) 

Ai 

The new physical position is r\ — biXi — bXo- In a numerical calculation the truncated b is exact to fioating point precision 
but contains an error because of the omitted orders; in a symbolic calculation b is known to order A™. 
The next step from ti — > t2 involves a similar update with b = b{ti — >■ ^2; Ai, 61) 

^ (E15) 

(Ef6) 

(E17) 

(l + '5i)(y)'-l (Ef8) 

<5„2 = ^-1. (Ef9) 

a20 

The new physical position is r2 = b2X2- This iterative scheme repeats for a total of TV steps. It ultimately yields an approxi- 
mation to the position at the final time denoted Nm ~ bMXN- 

A difference like (TV + l)m — TVm rnay be calculated numerically for various TV and m and the scaling fitted and inferred. 
In addition, one can approach the problem symbolically. To write TV^^ we need to expand the final result in powers of Aq. 
Note, for example, that Ai and 9i are known as expansions in Aq with coefficients that depend upon 9o. Perturbation-related 
quantities are re- written systematically in terms of initial quantities. For example, b{ti — > t2;Ai,6i) may be expanded in 
powers of Ai with coefficients depending upon 61. Next, all occurrences of Ai and 9i are replaced by expansions in powers 
of Ao and coefficients depending upon 9o. All terms up to and including A™ are retained in the final result. This procedure 
is systematically repeated until all quantities are expressed in terms of initial data. Finally, the difference (TV + f)™, — TV^ is 
calculated symbolically. Similar strategies allow construction of all the differences of interest. 

To make analytic progress assume that ft ~ P — 1 — {tf — to)/to « 1 is a small parameter. In a difference like 
(TV -I- l)m — TVm many "lower order" terms will coincide. Consider an ordering of terms by the powers of Ao (first) and by 
powers of ft (second). Define the leading-order difference to be the first non- vanishing term proportional to Ag/' for smallest 
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p and then smallest q. It is straightforward to apply this ordering to simplify the differences like (A^ + l)m — Nm- The leading 
order differences satisfy the following simple equalities 

I loo — Nm\ ~ IffiV.ml 
|iVm+l-iVml ~ IffiV.ml (E20) 

\{N+l)m — Nm\ ~ |3iV + l,m — gjV,m| 

where 



ffjv.m = KN,m COS sm 6* A /t 
iV- 
9 V 3 y V ^"+1 



^ — ^ ( ~'^\ { " 2 + m 



These differences can be compared with the numerical differences for which no expansion in ft is carried out. 

We verified the analytical scaling by the following numerical experiment. The parameters of the problem at the starting 
time to are Aq = 1/2, ^o = — 7r/4. The final time of interest tf is close to the initial time so that {tf — to) /to = 1/4. The m-th 
order Lagrangian approximation is evaluated at this fixed final time with successively increasing number of steps. Values of m 
from 1 to 4 and values of A'^ from 10 to 50 were considered. For geometric time steps (ti+i — ti)/ti = /3^'^ — 1 is independent 
of i and denoted 5t/t below. 

The results are plotted in figure |E1| The points indicate the numerical data points and the solid lines indicate the 



analytical functions defined in eq. (E20I. The numerical calculation was done with a high enough precision that even small 
errors of the order of 10~^^ are not contaminated by floating point errors. The agreement between the numerical experiment 
and the symbolic differences is very good. 

Thus, the scaling of the errors implies that for a small total time step, any finite order Lagrangian scheme will yield 
convergent results upon taking multiple steps. Conversely, for a fixed number of steps, a higher order Lagrangian calculation 
will give better results. 

It is useful to express the scaling in terms of the individual small step size 5t/t. Under the assumptions that [tf — to) /to 
is small, (f/ - io)/to ~ N5t/t. The scaling |loo - iV™| ~ iV"™A™+^((i/ - to)/to)"^+'^ can be re-written as |loo - N^\ ~ 
N{{tf ~ to)/to) ■ A'"+^((5f/f)'"+\ which can be interpreted as an error of {{tf - to)/to)A'"+^((5f/f)'"+^ per step. In the text, 
the quantity e — A5t/t is kept constant. For fixed initial and final times, the error scales oc A^e™"*"^. If A does not change 
appreciably then the error is oc Ae™. Convergence is attained when e — >■ 0. 
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Figure D3. Roots with rj = 2it plotted in the complex A plane for < 9 ^ w. These values of 9 correspond to those in figure [PT] The 
colour codes the complex phase of the roots (A = Ae**^). The real positive {(j> = 0) and negative {<j> = n) roots are shown in red and cyan 
respectively. The complex roots can have any colour other than these two and the bottom figure provides the coding. By comparison with 
figure [D1] one sees that all open models with real roots are cyan (negative); likewise all closed models with real roots are red (positive). 
Note, however, that there exist complex roots for both open and closed models. Coloured version online. 
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LogljN 
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Figure El. The three panels show the log of the errors |loo — Nm\, l-^Vm+i — Nm\ and \{N + l)m — Nm\ vs. A''. The final time tf is 
the same for all these comparisons. The dots correspond to the data generated by the numerical experiment and lines correspond to the 
analytical formulas given in cq. | |E20[ |. The lines from top to bottom correspond to m = 1, 2, 3, 4 respectively for the first and third panels 
and to m = 1, 2, 3 for the second panel. It is clear that for a fixed m, increasing the number of steps improves convergence. Conversely, 
for a fixed A'^, increasing the Lagrangian order m improves convergence. 
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